BEI and equity paper pre-analyses

1 Proposal

1.1 Original grant proposal for Census tract level analyses

1.1.1 Objective 1: Socio-spatial inequalities in urban interventions:

• To test if urban interventions tend to be located in low SES neighborhoods (H1a); we will run the following model: \[Log(𝔼(UI_{𝑖} \mid X_{𝑖})) = \alpha + \beta_{1} * SES_{𝑖}\] a Poisson regression model with \(UI\) = the number of Urban Interventions between 2016 and 2021, \(SES\) = the socio-economic status of the census track at 2016.
• To test the reduction in socio-economic inequalities in urban conditions (e.g. length and type of bicycle lanes, NDVI) between 2006 and 2020 (H1b), we will run the following model: \[𝔼(UC_{ij} \mid X_{ij}) = \beta_{0j} + \beta_{1j} ∗ SES + \beta_{2j} ∗ Time + \beta_{3j} ∗ SES ∗ Time + \epsilon_{ij}\] a multilevel linear model with \(UC\) = the Urban condition, \(SES\) = the socio-economic status of the census track.

1.1.2 Objective 2: Examining causal pathways and directionality of intervention/gentrification relation:

• To test if urban interventions occur before gentrification or if gentrification occurs before the implementation of new urban interventions, we will run two models : (1) \[Logit(𝔼(G_{i} \mid X_{i})) = \alpha + \beta_{1} ∗ UI_{i} + + \beta_{2} ∗ PG_{i} + \beta ∗ X_{control_{i}}\] a logistic regression model among low SES census tracts (2006) with \(G\) = Gentrification between 2006 and 2016, \(UI\) = the number of Urban Interventions between 2001 and 2006, and \(PG\) = SES status 2006; and (2) \[Logit(𝔼(UI_{i} \mid X_{i})) = \alpha + \beta_{1} ∗ G_{i} + \beta_{2} ∗ UC_{i} + \beta ∗ X_{control_{i}}\] a Poisson regression model, with \(UI\) = the number of Urban Interventions between 2016 and 2023, \(G\) = Gentrification between 2006 and 2016, and \(UC\) = Urban Conditions in 2016.

1.2 Paper objective

We focus on obj #1: are urban interventions tend to be located in low SES neighborhoods. Here, the urban interventions considered are bike lanes and canopy/tree coverage and low SES ~ high Pampalon deprivation index. In a second step, we will look at the variations of UI and SES and their association.

Data extraction and pre-analyses for the paper looking at BEI and equity. BEI comprise bike lanes and canopy changes while equity is measured through Pampalon deprivation index. (Partially adapted from original work on BEI bike lanes – see bike_lane_stats.R and ReadMe.md.)

Paper is available here.

General processing steps:

• Get CT 2016 boundaries (from Cancensus API)
• Get BEI interventions (canopy / bike lanes) for the selected years (2011 / 2015 / 2016 / 2019)
• Compute changes between years
  • for the original CT boundaries
  • for buffered CTs (250 / 500 / 750m)

In a second phase, these BEI changes will be linked to Pampalon index for 2016 (and 2011 ?)

1.3 Revision history

UPDATE 2021-12-02 Following discussion with @Yan, add normalized bike line changes:

• by CT/buffer area
• by street length within CT/buffer

UPDATE 2021-12-08 Following discussion with @Yan

• Keep only normalized variation of bike lane length
• Reverse relation: BEI metric = f(Equity metric (Pampalon / gentrification))
• Display association at baseline, then add delta over time
• Canopy : distinguish between high and low canopy

UPDATE 2022-01-14 Following discussion with @Ruben and @Yan

• Test model with UC 2011 when modeling UI
• Probably no pertinent to
• Try to define Year variable as continuous instead of category in LME
• Check if range of wSCOREMAT matches theoretical range of x-axis in graph (-.2 >> .2), might ponder the magnitude of the observed trend
• Test three-way interaction instead of 2 times two-way interactions (interactions might be complicated to interpret…)
• For gentrification and multilevel model, use GEE (logit)

2 Built Environment Intervention Extraction

2.1 Bike lane changes

We use data categorized by Philippe Apparicio’s team who manually identified bike lanes for each census year since 1991. For this study, we limit ourselves to 2016 and 2011 census years.

On top of the original CT boundaries, three levels of buffer have been applied to the CT – 250m, 500m & 750m. Then the same series of processing steps (see above) have been applied to the buffers.

# Bike lanes, from Ph. Apparicio
reseau <- st_read(dsn="data/ReseauCyclableFinal.gdb", layer = "Reseau") # Already in NAD83 / MTM zone 8

## Reading layer `Reseau' from data source 
##   `/Users/benoit/WORKSPACE/gentrification_BEI_equity/data/ReseauCyclableFinal.gdb' 
##   using driver `OpenFileGDB'
## Simple feature collection with 82166 features and 72 fields
## Geometry type: GEOMETRY
## Dimension:     XYZ, XYZM
## Bounding box:  xmin: 266985.5 ymin: 5029251 xmax: 320986.1 ymax: 5062652
## z_range:       zmin: 0 zmax: 43
## m_range:       mmin: 0 mmax: 43
## Projected CRS: NAD83 / MTM zone 8

bike_lane <- reseau %>%
  filter(An2016 == 1 | An2011 == 1) %>%
  select(IdRte, ClsRte, Zone, starts_with("An"), starts_with("Typo_")) %>%
  st_cast("MULTILINESTRING") # Get rid of a few MULTICURVE geometries

# CT boundaries for Montreal
CT16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='CT', geo_format = "sf") %>%
  filter(Type == "CT") %>%
  mutate(interact_aoi = CD_UID %in% c(2466, 2465, 2458)) %>% # Flag Montréal island, Laval and the South shore (Longueuil, St-Lambert, Brossard)
  st_transform(st_crs(bike_lane))

## Reading geo data from local cache.

CT11 <- get_census(dataset='CA11', regions=list(CMA='24462'), level='CT', geo_format = "sf") %>%
  filter(Type == "CT") %>%
  st_transform(st_crs(bike_lane))

## Reading geo data from local cache.

compute_bikelane_by_area <- function(sf_areas, year_fld, typo_fld) {
  # Compute length of bike lanes within each area.
  # --
  # Parameters:
  #   - sf_areas: sf class object defining the areas of interest, must have a GeoUID field
  #   - year_fld: field name specifying the year of interest, e.g. An2016
  #   - typo_fld: field name specifying the typology for the year of interest, e.g. Typo_2016
  
  year_fld <- enquo(year_fld)
  typo_fld <- enquo(typo_fld)
  
  # Compute intersection of bike lanes with areas
  bk <- bike_lane %>%
    filter(!!year_fld == 1) %>%
    st_intersection(sf_areas) %>%
    mutate(bike_lane_length = st_length(.)) %>%
    as.data.frame() %>%
    group_by(GeoUID, !!typo_fld) %>%
    summarise(bike_lane_length = sum(bike_lane_length)) %>%
    ungroup() %>%
    pivot_wider(names_from = !!typo_fld, names_prefix = "Bike_class", names_sort = TRUE,
                values_from = bike_lane_length, values_fill = units::set_units(0, m)) %>%
    mutate(Bike_lane_total = units::set_units(rowSums(select(., starts_with("Bike_class"))), m))
  
  # Merge back into original sf_areas
  bk <- sf_areas %>%
    left_join(bk)
  
  # Replace NA by 0, which occur in Bike_class length
  bk[is.na(bk)] <- 0
  
  return(bk)
}

compute_streetlength_by_area <- function(sf_areas) {
  # Compute length of streets within each area.
  # --
  # Parameters:
  #   - sf_areas: sf class object defining the areas of interest, must have a GeoUID field

  # Compute intersection of streets with areas
  bk <- reseau  %>%
    st_cast("MULTILINESTRING") %>% # Get rid of a few MULTICURVE geometries
    st_intersection(sf_areas) %>%
    mutate(street_length = st_length(.)) %>%
    as.data.frame() %>%
    group_by(GeoUID) %>%
    summarise(street_length = sum(street_length)) %>%
    ungroup()
  
  # Merge back into original sf_areas
  bk <- sf_areas %>%
    mutate(shape_area_km2 = units::set_units(st_area(.), 'km^2')) %>%
    left_join(bk)
  
  # Replace NA by 0
  bk[is.na(bk)] <- 0
  
  return(bk)
}


# Compute year 2016 and year 2011 bike lanes within 2016 CTs
# NB: contrary to the original work, we keep the same area of reference, i.e. 2016
bike_lane_by_CT16 <- compute_bikelane_by_area(CT16, An2016, Typo_2016)
bike_lane_by_CT11 <- compute_bikelane_by_area(CT16, An2011, Typo_2011)
bike_lane_by_CT06 <- compute_bikelane_by_area(CT16, An2006, Typo_2006)

# Compute buffers, with 3 radii and for each census year
radii <- c(250, 500, 750)

buf_CT16 <- lapply(radii, st_buffer, x=CT16)
names(buf_CT16) <- lapply(radii, function(b) {paste0("buf", b)})

# Compute bike length for each buffer/census year
buf_CT16_w_bike_length <- lapply(buf_CT16, compute_bikelane_by_area, year_fld=An2016, typo_fld=Typo_2016)
names(buf_CT16_w_bike_length) <- lapply(radii, function(b) {paste0("buf", b)})
buf_CT16_w_bike_length$original <- bike_lane_by_CT16

buf_CT11_w_bike_length <- lapply(buf_CT16, compute_bikelane_by_area, year_fld=An2011, typo_fld=Typo_2011)
names(buf_CT11_w_bike_length) <- lapply(radii, function(b) {paste0("buf", b)})
buf_CT11_w_bike_length$original <- bike_lane_by_CT11

buf_CT06_w_bike_length <- lapply(buf_CT16, compute_bikelane_by_area, year_fld=An2006, typo_fld=Typo_2006)
names(buf_CT06_w_bike_length) <- lapply(radii, function(b) {paste0("buf", b)})
buf_CT06_w_bike_length$original <- bike_lane_by_CT06

# Compute total street length within CT/buffer
street_length_by_CT16 <- compute_streetlength_by_area(CT16)
buf_CT16_street_length <- lapply(buf_CT16, compute_streetlength_by_area)
names(buf_CT16_street_length) <- lapply(radii, function(b) {paste0("buf", b)})
buf_CT16_street_length$original <- street_length_by_CT16

# Reorganize data to have all data in one dataframe
bike_lane_changes <- CT16 %>%
  left_join(select(as.data.frame(buf_CT16_street_length$original), GeoUID, street_length, shape_area_km2), by="GeoUID") %>%
  left_join(select(as.data.frame(buf_CT16_street_length$buf250), GeoUID, street_length, shape_area_km2), by="GeoUID", suffix=c(".ct", ".b250")) %>%
  left_join(select(as.data.frame(buf_CT16_street_length$buf500), GeoUID, street_length, shape_area_km2), by="GeoUID") %>%
  left_join(select(as.data.frame(buf_CT16_street_length$buf750), GeoUID, street_length, shape_area_km2), by="GeoUID", suffix=c(".b500", ".b750")) %>%
  left_join(select(as.data.frame(buf_CT16_w_bike_length$original), GeoUID, starts_with("Bike_")), by="GeoUID") %>%
  left_join(select(as.data.frame(buf_CT11_w_bike_length$original), GeoUID, starts_with("Bike_")), by="GeoUID", suffix=c(".2016ct", ".2011ct")) %>%
  left_join(select(as.data.frame(buf_CT16_w_bike_length$buf250), GeoUID, starts_with("Bike_")), by="GeoUID") %>%
  left_join(select(as.data.frame(buf_CT11_w_bike_length$buf250), GeoUID, starts_with("Bike_")), by="GeoUID", suffix=c(".2016b250", ".2011b250")) %>%
  left_join(select(as.data.frame(buf_CT16_w_bike_length$buf500), GeoUID, starts_with("Bike_")), by="GeoUID") %>%
  left_join(select(as.data.frame(buf_CT11_w_bike_length$buf500), GeoUID, starts_with("Bike_")), by="GeoUID", suffix=c(".2016b500", ".2011b500")) %>%
  left_join(select(as.data.frame(buf_CT16_w_bike_length$buf750), GeoUID, starts_with("Bike_")), by="GeoUID") %>%
  left_join(select(as.data.frame(buf_CT11_w_bike_length$buf750), GeoUID, starts_with("Bike_")), by="GeoUID", suffix=c(".2016b750", ".2011b750")) %>%
  left_join(select(as.data.frame(buf_CT06_w_bike_length$original), GeoUID, starts_with("Bike_")), by="GeoUID") %>%
  left_join(select(as.data.frame(buf_CT06_w_bike_length$buf250), GeoUID, starts_with("Bike_")), by="GeoUID", suffix=c(".2006ct", ".2006b250")) %>%
  left_join(select(as.data.frame(buf_CT06_w_bike_length$buf500), GeoUID, starts_with("Bike_")), by="GeoUID") %>%
  left_join(select(as.data.frame(buf_CT06_w_bike_length$buf750), GeoUID, starts_with("Bike_")), by="GeoUID", suffix=c(".2006b500", ".2006b750"))

# Compute ratio of bike lane vs street length
bike_lane_changes <- bike_lane_changes %>%
  mutate(Bike_lane.by.street.2011ct = 100 * Bike_lane_total.2011ct / street_length.ct,
         Bike_lane.by.street.2011b250 = 100 * Bike_lane_total.2011b250 / street_length.b250,
         Bike_lane.by.street.2011b500 = 100 * Bike_lane_total.2011b500 / street_length.b500,
         Bike_lane.by.street.2011b750 = 100 * Bike_lane_total.2011b750 / street_length.b750,
         Bike_lane.by.street.2016ct = 100 * Bike_lane_total.2016ct / street_length.ct,
         Bike_lane.by.street.2016b250 = 100 * Bike_lane_total.2016b250 / street_length.b250,
         Bike_lane.by.street.2016b500 = 100 * Bike_lane_total.2016b500 / street_length.b500,
         Bike_lane.by.street.2016b750 = 100 * Bike_lane_total.2016b750 / street_length.b750,
         Bike_lane.by.street.2006ct = 100 * Bike_lane_total.2006ct / street_length.ct,
         Bike_lane.by.street.2006b250 = 100 * Bike_lane_total.2006b250 / street_length.b250,
         Bike_lane.by.street.2006b500 = 100 * Bike_lane_total.2006b500 / street_length.b500,
         Bike_lane.by.street.2006b750 = 100 * Bike_lane_total.2006b750 / street_length.b750)

# Compute change between 2011 and 2016 (only for total bike lane length)
bike_lane_changes <- bike_lane_changes %>%
  mutate(Bike_lane_diff.2011.2016ct = Bike_lane_total.2016ct - Bike_lane_total.2011ct,
         Bike_lane_diff.2011.2016b250 = Bike_lane_total.2016b250 - Bike_lane_total.2011b250,
         Bike_lane_diff.2011.2016b500 = Bike_lane_total.2016b500 - Bike_lane_total.2011b500,
         Bike_lane_diff.2011.2016b750 = Bike_lane_total.2016b750 - Bike_lane_total.2011b750)

# Normalize bike lane change by (i) street length and (ii) area
bike_lane_changes <- bike_lane_changes %>%
  mutate(Bike_lane_diff.by.street.2011.2016ct = 100 * Bike_lane_diff.2011.2016ct / street_length.ct,
         Bike_lane_diff.by.street.2011.2016b250 = 100 * Bike_lane_diff.2011.2016b250 / street_length.b250,
         Bike_lane_diff.by.street.2011.2016b500 = 100 * Bike_lane_diff.2011.2016b500 / street_length.b500,
         Bike_lane_diff.by.street.2011.2016b750 = 100 * Bike_lane_diff.2011.2016b750 / street_length.b750,
         Bike_lane_diff.by.area.2011.2016ct = Bike_lane_diff.2011.2016ct / shape_area_km2.ct,
         Bike_lane_diff.by.area.2011.2016b250 = Bike_lane_diff.2011.2016b250 / shape_area_km2.b250,
         Bike_lane_diff.by.area.2011.2016b500 = Bike_lane_diff.2011.2016b500 / shape_area_km2.b500,
         Bike_lane_diff.by.area.2011.2016b750 = Bike_lane_diff.2011.2016b750 / shape_area_km2.b750)

# Save results
st_write(bike_lane_changes, dsn = "data/bike_length_changes.gpkg", delete_layer = TRUE)

## Deleting layer `bike_length_changes' using driver `GPKG'
## Writing layer `bike_length_changes' to data source 
##   `data/bike_length_changes.gpkg' using driver `GPKG'
## Writing 970 features with 140 fields and geometry type Multi Polygon.

# Clean up
rm(buf_CT16)
rm(bike_lane_by_CT11, bike_lane_by_CT16, bike_lane_by_CT06)
rm(buf_CT11_w_bike_length, buf_CT16_w_bike_length, buf_CT06_w_bike_length)
rm(street_length_by_CT16, buf_CT16_street_length, buf_CT06_street_length)

2.1.1 Illustration of bike lane metrics

Check output for one specific dataset (Census tracts 2016, no buffer)

2.1.1.1 Raw length

Bike lane length in 2016 within CTs, measured in meters

ggplot() +
   geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_total.2016ct)), lwd=0) +
  scale_fill_continuous(name = "Total length (m)")+ 
  labs(title = "Length of bike lanes within 2016 CTs", subtitle = "(INTERACT study area || for control only)")

2.1.1.2 Absolute length change

Absolute bike lane length change between 2011 and 2016, in meters

ggplot() +
   geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.2011.2016ct)), lwd=0) +
  scale_fill_gradient2(name = "Length (m)")+ 
  labs(title = "Changes in bike lane between 2011 and 2016", subtitle = "(INTERACT study area || CT level || for control only)")

2.1.1.3 Normalized length change (by street)

Relative bike lane length change between 2011 and 2016, normalized by street length within CT in 2016, expressed in %

\[Variation = \frac{Bike Lane_{2016}}{Street Length_{2016}} - \frac{Bike Lane_{2011}}{Street Length_{2016}}\]

ggplot() +
   geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.by.street.2011.2016ct)), lwd=0) +
  scale_fill_gradient2(name = "Variation (%)")+ 
  labs(title = "Changes in bike lane between 2011 and 2016, normalized by street", subtitle = "(INTERACT study area || CT level || for control only)")

2.1.1.4 Normalized length change (by area)

Relative bike lane length change between 2011 and 2016, normalized by CT area, expressed in \(\frac{km}{km^{2}}\)

\[Ratio = \frac{Bike Lane_{2016}}{CT Area_{2016}} - \frac{Bike Lane_{2011}}{CT Area_{2016}}\]

ggplot() +
   geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.by.area.2011.2016ct)), lwd=0) +
  scale_fill_gradient2(name = "Ratio (1/km)")+ 
  labs(title = "Changes in bike lane between 2011 and 2016, normalized by area", subtitle = "(INTERACT study area || CT level || for control only)")

2.1.1.5 Normalized length change (by street) NEW

After some discussion with Yan, we envision using a relative bike lane length ratio change between 2011 and 2016, normalized by the ratio in 2011, expressed in %

\[Variation = \frac{\frac{Bike Lane_{2016}}{Street Length_{2016}} - \frac{Bike Lane_{2011}}{Street Length_{2016}}}{\frac{Bike Lane_{2011}}{Street Length_{2016}}} = \frac{Bike Lane_{2016} - Bike Lane_{2011}}{Bike Lane_{2011}}\]

Problem: all CT with no bike lane in 2011 get a missing data, contrary to the original metric, which measured the absolute variation of the ratio of bike lane to street length.

.bike_lane_changes <- bike_lane_changes %>%
  mutate(Bike_lane_diff.by.street.relative.2011.2016ct = Bike_lane_diff.by.street.2011.2016ct / (Bike_lane_total.2011ct / street_length.ct))

ggplot() +
   geom_sf(data=filter(.bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.by.street.relative.2011.2016ct)), lwd=0) +
  scale_fill_gradient2(name = "Ratio")+ 
  labs(title = "Changes in bike lane ratio between 2011 and 2016, normalized by ratio in 2011", subtitle = "(INTERACT study area || CT level || for control only)")

2.1.2 Basic stats on each layer

2.1.2.1 Census Tracts

# CT level
ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.2011.2016ct)), binwidth = 250) + 
  xlab("Difference of bike lane length between 2016 and 2011 | CT level")

summary(bike_lane_changes$Bike_lane_diff.2011.2016ct)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1377.7     0.0     0.0   262.1   189.7 14463.8

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.street.2011.2016ct))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by street")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 252 rows containing non-finite values (stat_bin).

summary(bike_lane_changes$Bike_lane_diff.by.street.2011.2016ct)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -100.000    0.000    0.000    3.665    4.496  100.000      252

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.area.2011.2016ct))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by area")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(bike_lane_changes$Bike_lane_diff.by.area.2011.2016ct)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.1474  0.0000  0.0000  0.4237  0.1764  8.9852

ggplot(.bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.street.relative.2011.2016ct))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by ratio in 2011")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 461 rows containing non-finite values (stat_bin).

summary(.bike_lane_changes$Bike_lane_diff.by.street.relative.2011.2016ct)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
## -100.0000    0.0000    0.2199       Inf  102.7796       Inf       367

2.1.2.2 Buffers 250m

# buf250 level
ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.2011.2016b250)), binwidth = 250) + 
  xlab("Difference of bike lane length between 2016 and 2011 | buf 250m")

summary(bike_lane_changes$Bike_lane_diff.2011.2016b250)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -3452.5     0.0     0.0   646.7   933.8 20786.8

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.street.2011.2016b250))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by street")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 233 rows containing non-finite values (stat_bin).

summary(bike_lane_changes$Bike_lane_diff.by.street.2011.2016b250)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -38.5802   0.0000   0.8538   3.3168   4.9043  56.6588      233

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.area.2011.2016b250))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by area")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(bike_lane_changes$Bike_lane_diff.by.area.2011.2016b250)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.4622  0.0000  0.0000  0.4001  0.4306  6.2156

2.1.2.3 Buffers 500m

# buf500 level
ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.2011.2016b500)), binwidth = 250) + 
  xlab("Difference of bike lane length between 2016 and 2011 | buf 500m")

summary(bike_lane_changes$Bike_lane_diff.2011.2016b250)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -3452.5     0.0     0.0   646.7   933.8 20786.8

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.street.2011.2016b500))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by street")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 229 rows containing non-finite values (stat_bin).

summary(bike_lane_changes$Bike_lane_diff.by.street.2011.2016b500)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -26.891   0.000   1.385   3.298   5.323  22.121     229

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.area.2011.2016b500))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by area")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(bike_lane_changes$Bike_lane_diff.by.area.2011.2016b500)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.94397  0.00000  0.03393  0.38745  0.49370  4.00417

2.1.2.4 Buffers 750m

# buf750 level
ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.2011.2016b750)), binwidth = 250) + 
  xlab("Difference of bike lane length between 2016 and 2011 | buf 750m")

summary(bike_lane_changes$Bike_lane_diff.2011.2016b750)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -3697.5     0.0   578.6  1842.2  2928.6 29325.1

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.street.2011.2016b750))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by street")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 222 rows containing non-finite values (stat_bin).

summary(bike_lane_changes$Bike_lane_diff.by.street.2011.2016b750)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -4.7579  0.1418  1.7793  3.2961  5.2054 22.5448     222

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.area.2011.2016b750))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by area")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(bike_lane_changes$Bike_lane_diff.by.area.2011.2016b750)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.59760  0.00000  0.07793  0.38166  0.53774  3.95463

2.2 Canopy changes

Canopy changes is based on data produced by CMM, using multispectral aerial imagery and lidar. In order to sync the observations with the census years, we focus on 2011 and 2019 with one extra observation point in 2015.

The processing steps are similar to the ones for the bike lanes:

• Import the raster for each of the 3 years
• Compute proportion of canopy within each area level (CT and buffers) for the 3 years

# Codes du raster "espace vert"
# 0. No data (hors CMM)
# 1. NDVI < 0,3 et MNH < 3,0m = Minéral bas (route, stationnement, etc.)
# 2. NDVI < 0,3 et MNH ≥ 3,0m = Minéral haut (constructions)
# 3. NDVI ≥ 0,3 et MNH < 3,0m = Végétal bas (culture, gazon, etc.)
# 4. NDVI ≥ 0,3 et MNH ≥ 3,0m = Végétal haut (canopée)
# 5. Aquatique

# Load rasters into pg database for further processing
system('psql -d gentrif_bei -c "CREATE EXTENSION IF NOT EXISTS postgis"')
system('psql -d gentrif_bei -c "CREATE EXTENSION IF NOT EXISTS postgis_raster"')

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2019') IS NOT NULL;")) == 0) {
  system("raster2pgsql -s 32188 -I -C -M data/canopy/2019/*.tif -F -t 1000x1000 canopee2019 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2019' already imported") }

## PG Raster 'canopee2019' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2017') IS NOT NULL;")) == 0) {
  system("raster2pgsql -s 32188 -I -C -M data/canopy/2017/*.tif -F -t 1000x1000 canopee2017 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2017' already imported") }

## PG Raster 'canopee2017' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2015') IS NOT NULL;")) == 0) {
  system("raster2pgsql -s 32188 -I -C -M data/canopy/2015/*.tif -F -t 1000x1000 canopee2015 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2015' already imported") }

## PG Raster 'canopee2015' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2011') IS NOT NULL;")) == 0) {
  system("raster2pgsql -s 32188 -I -C -M data/canopy/2011/*.tif -F -t 1000x1000 canopee2011 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2011' already imported") }

## PG Raster 'canopee2011' already imported

# Resample to 10m as the original rasters have a 1m resolution, which is too high to allow for a swift processing
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2019_10m') IS NOT NULL;")) == 0) {
  system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2019 mode=2\" -r mode -tr 10 10 data/canopy/canopee2019_10m.tif")
  system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2019_10m.tif -F -t 100x100 canopee2019_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2019_10m' already imported") }

## PG Raster 'canopee2019_10m' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2017_10m') IS NOT NULL;")) == 0) {
  system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2017 mode=2\" -r mode -tr 10 10 data/canopy/canopee2017_10m.tif")
  system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2017_10m.tif -F -t 100x100 canopee2017_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2017_10m' already imported") }

## PG Raster 'canopee2017_10m' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2015_10m') IS NOT NULL;")) == 0) {
  system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2015 mode=2\" -r mode -tr 10 10 data/canopy/canopee2015_10m.tif")
  system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2015_10m.tif -F -t 100x100 canopee2015_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2015_10m' already imported") }

## PG Raster 'canopee2015_10m' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2011_10m') IS NOT NULL;")) == 0) {
  system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2011 mode=2\" -r mode -tr 10 10 data/canopy/canopee2011_10m.tif")
  system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2011_10m.tif -F -t 100x100 canopee2011_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2011_10m' already imported") }

## PG Raster 'canopee2011_10m' already imported

# Push CT16 to pg
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('ct16') IS NOT NULL;")) == 0) {
  CT16 %>%
    st_transform(crs = 32188) %>%
    st_write(con_bei, "ct16",
             layer_options = c("OVERWRITE=yes", "LAUNDER=true", "SPATIAL_INDEX=gist", "GEOMETRY_NAME=geom"))
  system("psql -d gentrif_bei -c 'CREATE INDEX ON  ct16 USING gist (geometry)'")
} else { message("PG Layer CT16 already imported") }

## PG Layer CT16 already imported

2.2.1 Extract % of green spaces at various scales

UPDATE 2021-12-08: following discussion with Yan, we decide to focus on extracting canopy years in sync with census years, i.e. 2011 and 2017, and 2015 as an intermediary year. (Previously, we were using the last available year, i.e. 2019), plus keep the option of looking high/low canopy separately

UPDATE 2022-02-09: adding the number of sq.m per inhabitant within census. when population of census equal or below 5 inhabitants.

2.2.1.1 Census Tracts

WITH cnt19 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2019_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee19 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2019_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2019
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2019
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2019
    FROM cnt19
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt17 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2017_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee17 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2017_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2017
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2017
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2017
    FROM cnt17
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt15 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2015_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee15 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2015_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2015
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2015
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2015
    FROM cnt15
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt11 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2011_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee11 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2011_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2011
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2011
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2011
    FROM cnt11
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
)
SELECT "GeoUID"
    ,st_area(geometry) ct_area_m2
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2011_by_pop, 0) END m2_esp_vert_2011_by_pop
    ,coalesce(pct_esp_vert_low_2011, 0) pct_esp_vert_low_2011
    ,coalesce(pct_esp_vert_high_2011, 0) pct_esp_vert_high_2011
    ,coalesce(pct_esp_vert_2011, 0) pct_esp_vert_2011
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2015_by_pop, 0) END m2_esp_vert_2015_by_pop
    ,coalesce(pct_esp_vert_low_2015, 0) pct_esp_vert_low_2015
    ,coalesce(pct_esp_vert_high_2015, 0) pct_esp_vert_high_2015
    ,coalesce(pct_esp_vert_2015, 0) pct_esp_vert_2015
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2017_by_pop, 0) END m2_esp_vert_2017_by_pop
    ,coalesce(pct_esp_vert_low_2017, 0) pct_esp_vert_low_2017
    ,coalesce(pct_esp_vert_high_2017, 0) pct_esp_vert_high_2017
    ,coalesce(pct_esp_vert_2017, 0) pct_esp_vert_2017
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2019_by_pop, 0) END m2_esp_vert_2019_by_pop
    ,coalesce(pct_esp_vert_low_2019, 0) pct_esp_vert_low_2019
    ,coalesce(pct_esp_vert_high_2019, 0) pct_esp_vert_high_2019
    ,coalesce(pct_esp_vert_2019, 0) pct_esp_vert_2019
FROM ct16
FULL JOIN canopee19 USING ("GeoUID")
FULL JOIN canopee17 USING ("GeoUID")
FULL JOIN canopee15 USING ("GeoUID")
FULL JOIN canopee11 USING ("GeoUID");

2.2.1.2 Buffers 250m

WITH ct16 AS (
    select "GeoUID", "Population", ST_Buffer(geometry, 250) geometry
    from ct16
),
cnt19 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2019_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee19 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2019_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2019
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2019
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2019
    FROM cnt19
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt17 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2017_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee17 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2017_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2017
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2017
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2017
    FROM cnt17
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt15 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2015_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee15 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2015_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2015
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2015
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2015
    FROM cnt15
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt11 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2011_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee11 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2011_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2011
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2011
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2011
    FROM cnt11
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
)
SELECT "GeoUID"
    ,st_area(geometry) ct_area_m2
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2011_by_pop, 0) END m2_esp_vert_2011_by_pop
    ,coalesce(pct_esp_vert_low_2011, 0) pct_esp_vert_low_2011
    ,coalesce(pct_esp_vert_high_2011, 0) pct_esp_vert_high_2011
    ,coalesce(pct_esp_vert_2011, 0) pct_esp_vert_2011
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2015_by_pop, 0) END m2_esp_vert_2015_by_pop
    ,coalesce(pct_esp_vert_low_2015, 0) pct_esp_vert_low_2015
    ,coalesce(pct_esp_vert_high_2015, 0) pct_esp_vert_high_2015
    ,coalesce(pct_esp_vert_2015, 0) pct_esp_vert_2015
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2017_by_pop, 0) END m2_esp_vert_2017_by_pop
    ,coalesce(pct_esp_vert_low_2017, 0) pct_esp_vert_low_2017
    ,coalesce(pct_esp_vert_high_2017, 0) pct_esp_vert_high_2017
    ,coalesce(pct_esp_vert_2017, 0) pct_esp_vert_2017
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2019_by_pop, 0) END m2_esp_vert_2019_by_pop
    ,coalesce(pct_esp_vert_low_2019, 0) pct_esp_vert_low_2019
    ,coalesce(pct_esp_vert_high_2019, 0) pct_esp_vert_high_2019
    ,coalesce(pct_esp_vert_2019, 0) pct_esp_vert_2019
FROM ct16
FULL JOIN canopee19 USING ("GeoUID")
FULL JOIN canopee17 USING ("GeoUID")
FULL JOIN canopee15 USING ("GeoUID")
FULL JOIN canopee11 USING ("GeoUID");

2.2.1.3 Buffers 500m

WITH ct16 AS (
    select "GeoUID", "Population", ST_Buffer(geometry, 500) geometry
    from ct16
),
cnt19 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2019_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee19 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2019_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2019
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2019
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2019
    FROM cnt19
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt17 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2017_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee17 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2017_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2017
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2017
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2017
    FROM cnt17
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt15 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2015_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee15 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2015_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2015
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2015
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2015
    FROM cnt15
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt11 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2011_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee11 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2011_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2011
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2011
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2011
    FROM cnt11
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
)
SELECT "GeoUID"
    ,st_area(geometry) ct_area_m2
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2011_by_pop, 0) END m2_esp_vert_2011_by_pop
    ,coalesce(pct_esp_vert_low_2011, 0) pct_esp_vert_low_2011
    ,coalesce(pct_esp_vert_high_2011, 0) pct_esp_vert_high_2011
    ,coalesce(pct_esp_vert_2011, 0) pct_esp_vert_2011
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2015_by_pop, 0) END m2_esp_vert_2015_by_pop
    ,coalesce(pct_esp_vert_low_2015, 0) pct_esp_vert_low_2015
    ,coalesce(pct_esp_vert_high_2015, 0) pct_esp_vert_high_2015
    ,coalesce(pct_esp_vert_2015, 0) pct_esp_vert_2015
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2017_by_pop, 0) END m2_esp_vert_2017_by_pop
    ,coalesce(pct_esp_vert_low_2017, 0) pct_esp_vert_low_2017
    ,coalesce(pct_esp_vert_high_2017, 0) pct_esp_vert_high_2017
    ,coalesce(pct_esp_vert_2017, 0) pct_esp_vert_2017
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2019_by_pop, 0) END m2_esp_vert_2019_by_pop
    ,coalesce(pct_esp_vert_low_2019, 0) pct_esp_vert_low_2019
    ,coalesce(pct_esp_vert_high_2019, 0) pct_esp_vert_high_2019
    ,coalesce(pct_esp_vert_2019, 0) pct_esp_vert_2019
FROM ct16
FULL JOIN canopee19 USING ("GeoUID")
FULL JOIN canopee17 USING ("GeoUID")
FULL JOIN canopee15 USING ("GeoUID")
FULL JOIN canopee11 USING ("GeoUID");

2.2.1.4 Buffers 750m

WITH ct16 AS (
    select "GeoUID", "Population", ST_Buffer(geometry, 750) geometry
    from ct16
),
cnt19 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2019_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee19 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2019_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2019
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2019
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2019
    FROM cnt19
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt17 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2017_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee17 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2017_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2017
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2017
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2017
    FROM cnt17
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt15 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2015_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee15 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2015_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2015
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2015
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2015
    FROM cnt15
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt11 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2011_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee11 AS (
    SELECT "GeoUID"
        ,round(10.*10. * sum(total) FILTER (WHERE value in (3, 4)) / NULLIF("Population", 0), 1) AS m2_esp_vert_2011_by_pop
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 2) AS pct_esp_vert_2011
        ,round(100. * sum(total) FILTER (WHERE value = 3) / sum(total), 2) AS pct_esp_vert_low_2011
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 2) AS pct_esp_vert_high_2011
    FROM cnt11
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
)
SELECT "GeoUID"
    ,st_area(geometry) ct_area_m2
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2011_by_pop, 0) END m2_esp_vert_2011_by_pop
    ,coalesce(pct_esp_vert_low_2011, 0) pct_esp_vert_low_2011
    ,coalesce(pct_esp_vert_high_2011, 0) pct_esp_vert_high_2011
    ,coalesce(pct_esp_vert_2011, 0) pct_esp_vert_2011
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2015_by_pop, 0) END m2_esp_vert_2015_by_pop
    ,coalesce(pct_esp_vert_low_2015, 0) pct_esp_vert_low_2015
    ,coalesce(pct_esp_vert_high_2015, 0) pct_esp_vert_high_2015
    ,coalesce(pct_esp_vert_2015, 0) pct_esp_vert_2015
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2017_by_pop, 0) END m2_esp_vert_2017_by_pop
    ,coalesce(pct_esp_vert_low_2017, 0) pct_esp_vert_low_2017
    ,coalesce(pct_esp_vert_high_2017, 0) pct_esp_vert_high_2017
    ,coalesce(pct_esp_vert_2017, 0) pct_esp_vert_2017
    ,CASE WHEN "Population" <= 5 THEN NULL ELSE COALESCE(m2_esp_vert_2019_by_pop, 0) END m2_esp_vert_2019_by_pop
    ,coalesce(pct_esp_vert_low_2019, 0) pct_esp_vert_low_2019
    ,coalesce(pct_esp_vert_high_2019, 0) pct_esp_vert_high_2019
    ,coalesce(pct_esp_vert_2019, 0) pct_esp_vert_2019
FROM ct16
FULL JOIN canopee19 USING ("GeoUID")
FULL JOIN canopee17 USING ("GeoUID")
FULL JOIN canopee15 USING ("GeoUID")
FULL JOIN canopee11 USING ("GeoUID");

2.3 Pampalon index

Get it here

pampalon <- read.xlsx("data/Canada2016Pampalon/A-MSDIData_Can2016_eng/1. EquivalenceTableCanada2016_ENG.xlsx", sheet = 2) %>%
  mutate(DA = as.character(DA)) %>%
  select(DA, SCOREMAT, SCORESOC)

# 2016 DA boundaries for Montreal
DA16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='DA', geo_format = "sf") %>%
  filter(Type == "DA") %>%
  st_transform(st_crs(bike_lane))

## Reading geo data from local cache.

pampalon <- DA16 %>%
  inner_join(pampalon, by = c("GeoUID" = "DA")) %>%
  as.data.frame()

# Get Pampalon 2006
pampalon06 <- read.xlsx("data/Canada2006Pampalon/A-MSDIData_Can2006_eng/1. CorrespondenceTable_Can2006_eng.xlsx", sheet = 2) %>%
  mutate(DA = as.character(DA)) %>%
  select(DA, DAPOP2006, SCOREMAT, SCORESOC)

# Get LUT DA2006 <-> DA2011 from StatCan
lut_da.1 <- read.csv("data/2011_92-156_DA_AD_txt/2011_92-156_DA_AD.txt", colClasses = "character", 
                     header = FALSE, col.names = c("DAUID2011.ADIDU2011", "DAUID2006.ADIDU2006", "DBUID2011", "DA_rel_flag")) %>%
  select(!c(DBUID2011, DA_rel_flag)) %>%
  unique()

# Link Pampalon 2011 to LUT and compute weighted mean of scores of Pampalon 2011
# NB: population numbers will diverge from  reality when more than one DA is merged into one DA of next census
pampalon06.11 <- pampalon06 %>%
  inner_join(lut_da.1, by = c("DA" = "DAUID2006.ADIDU2006")) %>%
  group_by(DAUID2011.ADIDU2011) %>%
  summarise(pop2006 = sum(DAPOP2006),
            SCOREMAT.06 = weighted.mean(SCOREMAT, DAPOP2006, na.rm = TRUE),
            SCORESOC.06 = weighted.mean(SCORESOC, DAPOP2006, na.rm = TRUE))

# Get Pampalon 2011
pampalon11 <- read.xlsx("data/Canada2011Pampalon/A-MSDIData_Can2011_eng/1. CorrespondenceTable_Can2011_eng.xlsx", sheet = 2) %>%
  mutate(DA = as.character(DA)) %>%
  select(DA, DAPOP2011, SCOREMAT, SCORESOC)

# Get LUT DA2011 <-> DA2016 from StatCan
lut_da <- read.csv("data/2016_92-156_DA_AD_csv/2016_92-156_DA_AD.csv", colClasses = "character") %>%
  select(!c(DBUID2016.IDIDU2016, DA_rel_flag.AD_ind_rel)) %>%
  unique()

# Link Pampalon 2011 to LUT, then to Pampalon 06 and finally compute weighted mean of scores of Pampalon 2011
pampalon11.16 <- pampalon11 %>%
  inner_join(lut_da, by = c("DA" = "DAUID2011.ADIDU2011")) %>%
  left_join(pampalon06.11, by =c("DA" = "DAUID2011.ADIDU2011")) %>%
  group_by(DAUID2016.ADIDU2016) %>%
  summarise(pop2011 = sum(DAPOP2011),
            SCOREMAT = weighted.mean(SCOREMAT, DAPOP2011, na.rm = TRUE),
            SCORESOC = weighted.mean(SCORESOC, DAPOP2011, na.rm = TRUE),
            SCOREMAT.06 = weighted.mean(SCOREMAT.06, pop2006, na.rm = TRUE),
            SCORESOC.06 = weighted.mean(SCORESOC.06, pop2006, na.rm = TRUE),
            pop2006 = sum(pop2006))

# Then link Pampalon 2011 to 2016
pampalon <- pampalon %>%
  left_join(pampalon11.16, by = c("GeoUID" = "DAUID2016.ADIDU2016"), suffix = c(".16", ".11"))

# Aggregate at the CT level
pampalon_CT <- pampalon %>%
  group_by(CT_UID) %>%
  summarise(wSCOREMAT.2016 = weighted.mean(SCOREMAT.16, Population, na.rm = TRUE),
            wSCORESOC.2016 = weighted.mean(SCORESOC.16, Population, na.rm = TRUE),
            wSCOREMAT.2011 = weighted.mean(SCOREMAT.11, pop2011, na.rm = TRUE),
            wSCORESOC.2011 = weighted.mean(SCORESOC.11, pop2011, na.rm = TRUE),
            wSCOREMAT.2006 = weighted.mean(SCOREMAT.06, pop2006, na.rm = TRUE),
            wSCORESOC.2006 = weighted.mean(SCORESOC.06, pop2006, na.rm = TRUE))

# Clean up
rm(lut_da, lut_da.1, pampalon11.16, pampalon06.11, pampalon11, pampalon06)

# Display map
.pampalon_CT_geom <- CT16 %>%
  left_join(pampalon_CT, by = c("GeoUID" = "CT_UID")) %>%
  filter(interact_aoi)

.pampalon_data <- bi_class(.pampalon_CT_geom, x = wSCOREMAT.2016, y = wSCORESOC.2016, style = "quantile", dim = 3)

## Warning in classInt::classIntervals(bins_x, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes

## Warning in classInt::classIntervals(bins_y, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes

.map <- ggplot() + 
  geom_sf(data = .pampalon_data, mapping = aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
  bi_scale_fill(pal = "DkBlue", dim = 3) +
  labs(title = "Pampalon: material and social deprivation index") + 
  theme(panel.background = element_rect(fill = "white"),
        #axis.ticks = element_blank(),
        #axis.text = element_blank(),
        panel.grid = element_line(color = "darkgray", size = 0.2))
.legend <- bi_legend(pal = "DkBlue",
                    dim = 3,
                    xlab = "Material ",
                    ylab = "Social ",
                    size = 8)
ggdraw() +
  draw_plot(.map, 0, 0, 1, 1) +
  draw_plot(.legend, 0.1, .7, 0.2, 0.2)

2.4 Gentrification

Using Ding metric computed on 5 year span.

# Load gentrified CTs, 5 year span (from repo gentrification_metrics)
ding <- list()
ding[["2016"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_16", quiet=TRUE) %>%
  filter(cma_uid_16 == "24462") %>%
  st_transform(st_crs(bike_lane))
ding[["2011"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_11", quiet=TRUE) %>%
  filter(cma_uid_11 == "24462") %>%
  st_transform(st_crs(bike_lane))
ding[["2006"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_06", quiet=TRUE) %>%
  filter(cma_uid_06 == "24462") %>%
  st_transform(st_crs(bike_lane))

.ding_map <- ding[["2016"]] %>%
  left_join(select(as.data.frame(CT16), GeoUID, interact_aoi), by = c("ct_uid_16" = "GeoUID")) %>%
  filter(interact_aoi)

ggplot(data = .ding_map) + 
  geom_sf(aes(fill = gentrified_2016_2011, colour=gentrifiable_2011)) +
  scale_fill_manual(values = c("gray", "red", "darkgray"), name = "Gentrified in 2016") +
  scale_colour_manual(values = c("darkgray", "darkred", "darkgray"), name = "Gentrifiable in 2011") +
  labs(title = "Census tract gentrification status in 2016")

2.5 Equity metrics | low-income population and visible minority

Introduced here as a proposition, nothing acted (2022-02-04)

# Visible Minority
# - v_CA16_3954: Total - Visible minority for the population in private households - 25% sample data (Total)
# - v_CA16_3957: Total visible minority population (Total)

# Low income (LIM-AT)
# - v_CA16_2540: Prevalence of low income based on the Low-income measure, after tax (LIM-AT) (%) (Total)
equity_ct16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='CT', geo_format = "sf",
                          vectors = c("v_CA16_3954", "v_CA16_3957", "v_CA16_2540")) %>%
  filter(Type == "CT") %>%
  transmute(CT_UID = GeoUID,
            vis_minority_2016 = `v_CA16_3957: Total visible minority population` / `v_CA16_3954: Total - Visible minority for the population in private households - 25% sample data` * 100,
            low_income_2016 = `v_CA16_2540: Prevalence of low income based on the Low-income measure, after tax (LIM-AT) (%)`)

## Reading vectors data from local cache.

## Reading geo data from local cache.

# Visible Minority
# - v_CA11N_457: CA 2011 NHS, Total population in private households by visible minority (Total)
# - v_CA11N_460: CA 2011 NHS, Total population in private households by visible minority, Total visible minority population (Total)

# Low income (LIM-AT)
# - v_CA11N_2606: CA 2011 NHS, Prevalence of low income in 2010 based on after-tax low-income measure % (Total)
equity_ct11 <- get_census(dataset='CA11', regions=list(CMA='24462'), level='CT', geo_format = "sf",
                          vectors = c("v_CA11N_457", "v_CA11N_460", "v_CA11N_2606")) %>%
  filter(Type == "CT") %>%
  transmute(CT_UID = GeoUID,
            vis_minority_2011 = `v_CA11N_460: Total visible minority population` / `v_CA11N_457: Total population in private households by visible minority` * 100,
            low_income_2011 = `v_CA11N_2606: Prevalence of low income in 2010 based on after-tax low-income measure %`)

## Reading vectors data from local cache.
## Reading geo data from local cache.

# Visible Minority
# - v_CA06_1302: Total population by visible minority groups
# - v_CA06_1303: Total population by visible minority groups, Total visible minority population

# Low income (LIM-AT)
# - v_TX2006_551: After-tax low income status of tax filers and dependents (census family low income measure, CFLIM-AT) for couple and lone parent families by family composition, 2006 | All family units | Persons in Low Income | % - Total
equity_ct06 <- get_census(dataset='CA06', regions=list(CMA='24462'), level='CT', geo_format = "sf",
                          vectors = c("v_CA06_1302", "v_CA06_1303", "v_TX2006_551")) %>%
  filter(Type == "CT") %>%
  transmute(CT_UID = GeoUID,
            vis_minority_2006 = `v_CA06_1303: Total visible minority population` / `v_CA06_1302: Total population by visible minority groups - 20% sample data` * 100,
            low_income_2006 = `v_TX2006_551: % - Total`)

## Reading vectors data from local cache.
## Reading geo data from local cache.

equity_ct <- st_join(equity_ct16, equity_ct11, left=TRUE, largest=TRUE, suffix=c("", "_2011")) %>% # join on largest overlap, to overcome mismatch in CT UID
  st_join(equity_ct06, left=TRUE, largest=TRUE, suffix=c("", "_2006")) %>%
  data.frame()

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

# cleanup
rm(equity_ct11, equity_ct16, equity_ct06)

# Display map
.equity_CT_geom <- CT16 %>%
  left_join(equity_ct, by = c("GeoUID" = "CT_UID")) %>%
  filter(interact_aoi)

.equity_data <- bi_class(.equity_CT_geom, x = vis_minority_2016, y = low_income_2016, style = "quantile", dim = 3)

## Warning in classInt::classIntervals(bins_x, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes

## Warning in classInt::classIntervals(bins_y, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes

.map <- ggplot() + 
  geom_sf(data = .equity_data, mapping = aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
  bi_scale_fill(pal = "Brown", dim = 3) +
  labs(title = "Equity metrics: % of visible minority and % of low-income household") + 
  theme(panel.background = element_rect(fill = "white"),
        #axis.ticks = element_blank(),
        #axis.text = element_blank(),
        panel.grid = element_line(color = "darkgray", size = 0.2))
.legend <- bi_legend(pal = "Brown",
                    dim = 3,
                    xlab = "Vis. Minority ",
                    ylab = "Low-Income ",
                    size = 8)
ggdraw() +
  draw_plot(.map, 0, 0, 1, 1) +
  draw_plot(.legend, 0.1, .7, 0.2, 0.2)

2.6 Build complete dataset

All variables + outcome linked at the CT level

.bike_lane_changes <- bike_lane_changes %>%
  as.data.frame() %>%
  select(GeoUID, ends_with("ct", ignore.case = FALSE), ends_with("b250", ignore.case = FALSE), ends_with("b500", ignore.case = FALSE), ends_with("b750", ignore.case = FALSE)) %>%
  select(GeoUID, starts_with("Bike_lane")) # Drop individual category lane length

bei_df <- CT16 %>%
  as.data.frame() %>%
  transmute(CT_UID = GeoUID,
            CD_UID = CD_UID,
            CSD_UID = CSD_UID,
            interact_aoi = interact_aoi,
            Population = Population) %>%
  left_join(pampalon_CT, by="CT_UID") %>%
  left_join(select(as.data.frame(ding$`2016`), ct_uid_16, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_16")) %>%
  left_join(select(as.data.frame(ding$`2011`), ct_uid_11, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_11")) %>%
  left_join(select(as.data.frame(ding$`2006`), ct_uid_06, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_06")) %>%
  left_join(select(as.data.frame(equity_ct), !c("geometry", "CT_UID_2011", "CT_UID_2006")), by="CT_UID") %>%
  left_join(.bike_lane_changes, by=c("CT_UID" = "GeoUID")) %>%
  left_join(as.data.frame(esp_vert_ct), by=c("CT_UID" = "GeoUID")) %>%
  left_join(as.data.frame(esp_vert_buf250), by=c("CT_UID" = "GeoUID"), suffix=c("ct", "b250")) %>%
  left_join(as.data.frame(esp_vert_buf500), by=c("CT_UID" = "GeoUID")) %>%
  left_join(as.data.frame(esp_vert_buf750), by=c("CT_UID" = "GeoUID"), suffix=c("b500", "b750"))
  
  
head(bei_df)

write.csv(bei_df, "data/_results/bei_equity.csv", na="", row.names = FALSE)

Included variables:

• Census Tracts variables
  • CT_UID: 2016 Census Tract ID
  • CD_UID: 2016 Census Division
  • CSD_UID: 2016 Census Subdivision
  • interact_aoi: Does CT belong to INTERACT study area?
  • Population: 2016 Population within CT
  • ct_area_m2.{ct|b{250|500|750}}: Area of CT or buffer of 250, 500 or 750m radius around CT, in square meters
• Gentrification metrics
  • gentrified_2016_2011: Is the CT gentrified in 2016?
  • gentrifiable_2011: Is the CT candidate to gentrification in 2011?
  • gentrified_2011_2006: Is the CT gentrified in 2011
  • gentrifiable_2006: Is the CT candidate to gentrification in 2006
  • gentrified_2006_2001: Is the CT gentrified in 2006
  • gentrifiable_2001: Is the CT candidate to gentrification in 2001
• Pampalon’s metrics
  • wSCOREMAT.2016: Social deprivation index in 2016 (population weighted)
  • wSCORESOC.2016: Material deprivation index in 2016 (population weighted)
  • wSCOREMAT.2011: Social deprivation index in 2011 (population weighted)
  • wSCORESOC.2011: Material deprivation index in 2011 (population weighted)
  • wSCOREMAT.2006: Social deprivation index in 2006 (population weighted)
  • wSCORESOC.2006: Material deprivation index in 2006 (population weighted)
• Social profile
  • vis_minority_2016: % of visible minority in CT 2016
  • low_income_2016: prevalence of low income in CT 2016
  • vis_minority_2011: % of visible minority in CT 2011
  • low_income_2011: prevalence of low income in CT 2011
  • vis_minority_2006: % of visible minority in CT 2006
  • low_income_2006: prevalence of low income in CT 2006
• Bike lane density
  • Bike_lane_total.{2016|2011}{ct|b{250|500|750}}: total length of bike lanes, in 2016 or 2011, within CT or buffer of 250, 500 or 750m radius
  • Bike_lane.by.street.{2016|2011}{ct|b{250|500|750}}: % of bike lanes compared to streets, in 2016 or 2011, within CT or buffer of 250, 500 or 750m radius
  • Bike_lane_diff.2011.2016{ct|b{250|500|750}}: change in total length of bike lanes between 2011 and 2011, within CT or buffer of 250, 500 or 750m radius
  • Bike_lane_diff.by.street.2011.2016{ct|b{250|500|750}}: change in total length of bike lanes between 2011 and 2011, normalized by street length, within CT or buffer of 250, 500 or 750m radius
  • Bike_lane_diff.by.area.2011.2016{ct|b{250|500|750}}: change in total length of bike lanes between 2011 and 2011, normalized by area, within CT or buffer of 250, 500 or 750m radius
• Green spaces
  • pct_esp_vert_{2011|2015|2019}.{ct|b{250|500|750}}: % of green space in 2011, 2015 or 2019 within CT or buffer of 250, 500 or 750m radius
  • pct_esp_vert_{low|high}_{2011|2015|2019}.{ct|b{250|500|750}}: same as above, except for grass (low) and tree (high)
  • pct_esp_vert_diff{2011|2015}.{2015|2019}.{ct|b{250|500|750}}: change in % of green space between 2011 and 2015, 2011 and 2019 as well as 2011 and 2019, within CT or buffer of 250, 500 or 750m radius

3 Preliminary analyses

3.1 INTERACT study area

INTERACT study area ~ Montréal, Laval, Longueuil, Brossard, St-Lambert

3.1.1 SES variable distribution

.bei_df_long <- bei_df %>% 
  filter(interact_aoi) %>%
  units::drop_units() %>% 
  select(CT_UID, CD_UID, starts_with("wSCORE")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name) #, scales = "free")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 86 rows containing non-finite values (stat_bin).

.bei_df_long <- bei_df %>% 
  filter(interact_aoi) %>%
  units::drop_units() %>% 
  select(CT_UID, CD_UID, starts_with("vis_minority"), starts_with("low_income")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name) #, scales = "free")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 85 rows containing non-finite values (stat_bin).

.bei_df_long <- bei_df %>% 
  filter(interact_aoi) %>%
  units::drop_units() %>% 
  select(CT_UID, CD_UID, starts_with("gentrif")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_bar() + 
  facet_wrap(~name) #, scales = "free", ncol = 3)

# Testing correlation between numeric variables
corrplot::corrplot(cor(select(bei_df, starts_with("wSCORE"), starts_with("vis_"), starts_with("low_")), use = "complete.obs"), type="upper", order="hclust")

# ANOVA for gentrified variables (see https://stats.stackexchange.com/questions/119835/correlation-between-a-nominal-iv-and-a-continuous-dv-variable)
heplots::etasq(aov(gentrified_2016_2011 ~ wSCOREMAT.2011, data = bei_df))

heplots::etasq(aov(gentrified_2016_2011 ~ low_income_2011, data = bei_df))

heplots::etasq(aov(gentrified_2016_2011 ~ vis_minority_2011, data = bei_df))

3.1.2 BEI variable distributions

3.1.2.1 Census Tracts

.bei_df_long <- bei_df %>% 
  filter(interact_aoi) %>%
  units::drop_units() %>% 
  select(CT_UID, CD_UID, matches("^Bike_lane.*ct$"), matches("^pct_esp_vert_diff.*ct$")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name, scales = "free", ncol = 3)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

3.1.2.2 Buffers 250m

.bei_df_long <- bei_df %>% 
  filter(interact_aoi) %>%
  units::drop_units() %>% 
  select(CT_UID, CD_UID, matches("^Bike_lane.*b250$"), matches("^pct_esp_vert_diff.*b250$")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name, scales = "free", ncol = 3)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

3.1.2.3 Buffers 500m

.bei_df_long <- bei_df %>% 
  filter(interact_aoi) %>%
  units::drop_units() %>% 
  select(CT_UID, CD_UID, matches("^Bike_lane.*b500$"), matches("^pct_esp_vert_diff.*b500$")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name, scales = "free", ncol = 3)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

3.1.2.4 Buffers 750m

.bei_df_long <- bei_df %>% 
  filter(interact_aoi) %>%
  units::drop_units() %>% 
  select(CT_UID, CD_UID, matches("^Bike_lane.*b750$"), matches("^pct_esp_vert_diff.*b750$")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name, scales = "free", ncol = 3)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

4 Association between SES and Urban Conditions at baseline (2011)

Looking at objective #1 | do urban interventions tend to be located in low SES neighborhoods?. We look at \[Urban Condition_{2011} = f(SES_{2011})\] as well as \[Urban Condition_{2011} = f(Gentrification_{2011 \to 2016})\]

Here \(UrbanCondition\) means the state of the urban environment features, such as length of bike lanes, greenness coverage, etc. This needs to be distinguished from \(UrbanIntervention\), which accounts for the changes in the \(UrbanConditions\) (see below).

4.1 INTERACT study area

# keep only interact CT
bei_df_aoi <- filter(bei_df, interact_aoi)

4.1.1 UC vs Pampalon | material

4.1.1.1 Bike lane length

Bike lane ratio to streets (in %)

4.1.1.1.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane.by.street.2011ct, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane.by.street.2011ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.200  -7.709  -2.064   5.143  59.255 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      8.8868     0.3611  24.613  < 2e-16 ***
## wSCOREMAT.2011 -27.1407     8.8021  -3.083  0.00213 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.409 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01363,    Adjusted R-squared:  0.0122 
## F-statistic: 9.508 on 1 and 688 DF,  p-value: 0.002128

4.1.1.1.2 Buffers 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane.by.street.2011b250, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane.by.street.2011b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.4394  -5.3193  -0.8782   3.8947  24.9430 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      8.9081     0.2544  35.019  < 2e-16 ***
## wSCOREMAT.2011 -33.4497     6.2014  -5.394 9.49e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.629 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04057,    Adjusted R-squared:  0.03918 
## F-statistic: 29.09 on 1 and 688 DF,  p-value: 9.492e-08

4.1.1.1.3 Buffers 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane.by.street.2011b500, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane.by.street.2011b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.7672  -3.9200  -0.7691   3.1333  20.7230 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      8.9691     0.2162  41.488  < 2e-16 ***
## wSCOREMAT.2011 -31.0529     5.2703  -5.892 5.97e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.634 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04804,    Adjusted R-squared:  0.04665 
## F-statistic: 34.72 on 1 and 688 DF,  p-value: 5.971e-09

4.1.1.1.4 Buffers 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane.by.street.2011b750, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane.by.street.2011b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5371  -3.5685  -0.5231   2.4075  18.3871 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      8.9986     0.1938  46.431  < 2e-16 ***
## wSCOREMAT.2011 -26.8673     4.7248  -5.686 1.92e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.05 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04489,    Adjusted R-squared:  0.0435 
## F-statistic: 32.34 on 1 and 688 DF,  p-value: 1.918e-08

4.1.1.2 Canopy

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %)

4.1.1.2.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_2011ct, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_2011ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.859 -11.624  -0.331   8.977  47.657 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      37.9617     0.5933  63.988  < 2e-16 ***
## wSCOREMAT.2011 -109.4130    14.4631  -7.565 1.25e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.46 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.07679,    Adjusted R-squared:  0.07545 
## F-statistic: 57.23 on 1 and 688 DF,  p-value: 1.247e-13

4.1.1.2.2 Buffer 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_2011b250, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_2011b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.973  -8.235  -0.954   7.170  41.284 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     36.0871     0.4897  73.690  < 2e-16 ***
## wSCOREMAT.2011 -90.8096    11.9387  -7.606 9.29e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.76 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.07757,    Adjusted R-squared:  0.07623 
## F-statistic: 57.86 on 1 and 688 DF,  p-value: 9.289e-14

4.1.1.2.3 Buffer 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_2011b500, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_2011b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.955  -8.076  -0.976   7.123  40.579 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     36.1464     0.4645  77.820  < 2e-16 ***
## wSCOREMAT.2011 -83.2220    11.3236  -7.349 5.66e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.1 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.07279,    Adjusted R-squared:  0.07145 
## F-statistic: 54.01 on 1 and 688 DF,  p-value: 5.662e-13

4.1.1.2.4 Buffer 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_2011b750, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_2011b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.621  -8.282  -0.894   7.305  39.867 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      36.211      0.451  80.287  < 2e-16 ***
## wSCOREMAT.2011  -76.126     10.995  -6.923 1.01e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.75 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06513,    Adjusted R-squared:  0.06377 
## F-statistic: 47.93 on 1 and 688 DF,  p-value: 1.013e-11

4.1.1.3 Canopy (trees)

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %)

4.1.1.3.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_high_2011ct, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_high_2011ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.954  -5.492  -0.802   4.630  44.704 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     18.8069     0.3379   55.66   <2e-16 ***
## wSCOREMAT.2011 -96.7463     8.2370  -11.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.805 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.167,  Adjusted R-squared:  0.1658 
## F-statistic:   138 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.1.3.2 Buffer 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_high_2011b250, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_high_2011b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.9166  -4.2910  -0.8451   3.4739  27.8788 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     17.6539     0.2668   66.17   <2e-16 ***
## wSCOREMAT.2011 -81.7658     6.5045  -12.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.953 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1868, Adjusted R-squared:  0.1856 
## F-statistic:   158 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.1.3.3 Buffer 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_high_2011b500, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_high_2011b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.4603  -4.0485  -0.9957   3.5770  26.5204 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      17.738      0.250   70.96   <2e-16 ***
## wSCOREMAT.2011  -76.880      6.094  -12.62   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.515 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1879, Adjusted R-squared:  0.1867 
## F-statistic: 159.1 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.1.3.4 Buffer 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_high_2011b750, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_high_2011b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.6658  -3.9100  -0.6752   3.4947  25.5210 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     17.7572     0.2377   74.69   <2e-16 ***
## wSCOREMAT.2011 -72.9660     5.7960  -12.59   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.196 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1872, Adjusted R-squared:  0.186 
## F-statistic: 158.5 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.2 UC vs Visible minority

4.1.2.1 Bike lane length

Bike lane ratio to streets (in %)

4.1.2.1.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane.by.street.2011ct, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane.by.street.2011ct ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011ct ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.713  -7.452  -2.008   4.863  59.699 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       11.07967    0.66409  16.684  < 2e-16 ***
## vis_minority_2011 -0.09265    0.02192  -4.226  2.7e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.289 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0253, Adjusted R-squared:  0.02388 
## F-statistic: 17.86 on 1 and 688 DF,  p-value: 2.702e-05

4.1.2.1.2 Buffers 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane.by.street.2011b250, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane.by.street.2011b250 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b250 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.1067  -4.7893  -0.6668   3.8226  22.6792 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       11.54371    0.46311  24.926  < 2e-16 ***
## vis_minority_2011 -0.11042    0.01529  -7.222 1.36e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.478 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.07046,    Adjusted R-squared:  0.06911 
## F-statistic: 52.15 on 1 and 688 DF,  p-value: 1.365e-12

4.1.2.1.3 Buffers 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane.by.street.2011b500, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane.by.street.2011b500 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b500 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.6525  -3.7275  -0.6147   3.0419  20.5437 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       11.50117    0.39113  29.405  < 2e-16 ***
## vis_minority_2011 -0.10577    0.01291  -8.191 1.26e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.471 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.08885,    Adjusted R-squared:  0.08752 
## F-statistic: 67.09 on 1 and 688 DF,  p-value: 1.264e-15

4.1.2.1.4 Buffers 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane.by.street.2011b750, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane.by.street.2011b750 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b750 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.6227  -3.3596  -0.4224   2.5213  17.7443 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       11.43262    0.34878  32.779   <2e-16 ***
## vis_minority_2011 -0.10094    0.01152  -8.766   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.879 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1005, Adjusted R-squared:  0.09915 
## F-statistic: 76.83 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.2.2 Canopy

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %)

4.1.2.2.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_2011ct, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_2011ct ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011ct ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.420 -12.261  -0.970   9.559  50.400 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       42.98501    1.12478  38.216  < 2e-16 ***
## vis_minority_2011 -0.22010    0.03713  -5.927 4.87e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.73 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04858,    Adjusted R-squared:  0.0472 
## F-statistic: 35.13 on 1 and 688 DF,  p-value: 4.873e-09

4.1.2.2.2 Buffer 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_2011b250, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_2011b250 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b250 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.184  -9.310  -1.106   7.525  44.699 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       39.79001    0.93414  42.595  < 2e-16 ***
## vis_minority_2011 -0.16391    0.03084  -5.315 1.45e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.07 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.03944,    Adjusted R-squared:  0.03804 
## F-statistic: 28.25 on 1 and 688 DF,  p-value: 1.445e-07

4.1.2.2.3 Buffer 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_2011b500, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_2011b500 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b500 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.877  -9.260  -1.111   7.543  41.841 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        39.2416     0.8876  44.212  < 2e-16 ***
## vis_minority_2011  -0.1386     0.0293  -4.732 2.71e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.42 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.03151,    Adjusted R-squared:  0.03011 
## F-statistic: 22.39 on 1 and 688 DF,  p-value: 2.706e-06

4.1.2.2.4 Buffer 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_2011b750, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_2011b750 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b750 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.932  -9.202  -0.917   7.572  40.308 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       38.87252    0.86145  45.124  < 2e-16 ***
## vis_minority_2011 -0.12045    0.02844  -4.235  2.6e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.05 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.02541,    Adjusted R-squared:  0.02399 
## F-statistic: 17.93 on 1 and 688 DF,  p-value: 2.597e-05

4.1.2.3 Canopy (trees)

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %)

4.1.2.3.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_high_2011ct, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_high_2011ct ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011ct ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.484  -6.304  -1.788   4.737  49.585 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       21.84613    0.67054  32.580  < 2e-16 ***
## vis_minority_2011 -0.13945    0.02214  -6.299 5.34e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.379 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05453,    Adjusted R-squared:  0.05315 
## F-statistic: 39.68 on 1 and 688 DF,  p-value: 5.342e-10

4.1.2.3.2 Buffer 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_high_2011b250, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_high_2011b250 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b250 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.891  -5.227  -1.552   3.870  30.072 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       19.68093    0.54102  36.377   <2e-16 ***
## vis_minority_2011 -0.09617    0.01786  -5.384    1e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.568 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04043,    Adjusted R-squared:  0.03903 
## F-statistic: 28.99 on 1 and 688 DF,  p-value: 1.001e-07

4.1.2.3.3 Buffer 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_high_2011b500, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_high_2011b500 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b500 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.866  -4.947  -1.760   3.788  30.662 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       19.21085    0.51074  37.614  < 2e-16 ***
## vis_minority_2011 -0.07341    0.01686  -4.353 1.54e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.144 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.02681,    Adjusted R-squared:  0.02539 
## F-statistic: 18.95 on 1 and 688 DF,  p-value: 1.544e-05

4.1.2.3.4 Buffer 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_high_2011b750, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_high_2011b750 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b750 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.188  -4.538  -1.596   3.704  28.055 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       18.90064    0.48742  38.777  < 2e-16 ***
## vis_minority_2011 -0.05982    0.01609  -3.718 0.000218 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.818 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01969,    Adjusted R-squared:  0.01827 
## F-statistic: 13.82 on 1 and 688 DF,  p-value: 0.0002175

4.1.3 UC vs gentrified CT

Gentrified CT between 2011 and 2016

4.1.3.1 Bike lane length

Bike lane ratio to streets (in %)

4.1.3.1.1 Census tract level

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=Bike_lane.by.street.2011ct, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(Bike_lane.by.street.2011ct, na.rm = TRUE),
    sd = sd(Bike_lane.by.street.2011ct, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(Bike_lane.by.street.2011ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)
## gentrified_2016_2011   1      5    4.95   0.056  0.813
## Residuals            688  60902   88.52               
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(Bike_lane.by.street.2011ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011ct ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.760 -8.575 -1.815  5.186 55.344 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                8.7600     0.4285  20.441   <2e-16 ***
## gentrified_2016_2011TRUE  -0.1846     0.7805  -0.236    0.813    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.409 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  8.126e-05,  Adjusted R-squared:  -0.001372 
## F-statistic: 0.05591 on 1 and 688 DF,  p-value: 0.8131

4.1.3.1.2 Buffers 250m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=Bike_lane.by.street.2011b250, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(Bike_lane.by.street.2011b250, na.rm = TRUE),
    sd = sd(Bike_lane.by.street.2011b250, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(Bike_lane.by.street.2011b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)
## gentrified_2016_2011   1      2    1.67   0.037  0.847
## Residuals            688  31057   45.14               
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(Bike_lane.by.street.2011b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b250 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7878 -4.9870 -0.6935  4.1027 23.9655 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                8.6805     0.3060  28.365   <2e-16 ***
## gentrified_2016_2011TRUE   0.1073     0.5574   0.192    0.847    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.719 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  5.382e-05,  Adjusted R-squared:  -0.0014 
## F-statistic: 0.03703 on 1 and 688 DF,  p-value: 0.8475

4.1.3.1.3 Buffers 500m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=Bike_lane.by.street.2011b500, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(Bike_lane.by.street.2011b500, na.rm = TRUE),
    sd = sd(Bike_lane.by.street.2011b500, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(Bike_lane.by.street.2011b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)
## gentrified_2016_2011   1     22   21.78   0.664  0.416
## Residuals            688  22580   32.82               
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(Bike_lane.by.street.2011b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b500 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.060 -4.157 -0.576  3.120 19.585 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                8.6728     0.2609  33.236   <2e-16 ***
## gentrified_2016_2011TRUE   0.3871     0.4753   0.815    0.416    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.729 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0009635,  Adjusted R-squared:  -0.0004886 
## F-statistic: 0.6635 on 1 and 688 DF,  p-value: 0.4156

4.1.3.1.4 Buffers 750m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=Bike_lane.by.street.2011b750, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(Bike_lane.by.street.2011b750, na.rm = TRUE),
    sd = sd(Bike_lane.by.street.2011b750, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(Bike_lane.by.street.2011b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)
## gentrified_2016_2011   1     41   40.84   1.547  0.214
## Residuals            688  18163   26.40               
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(Bike_lane.by.street.2011b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane.by.street.2011b750 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.2152 -3.8001 -0.7224  2.6198 19.1075 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                8.6850     0.2340  37.110   <2e-16 ***
## gentrified_2016_2011TRUE   0.5302     0.4263   1.244    0.214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.138 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.002243,   Adjusted R-squared:  0.0007931 
## F-statistic: 1.547 on 1 and 688 DF,  p-value: 0.214

4.1.3.2 Canopy

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %)

4.1.3.2.1 Census tract level

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_2011ct, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_2011ct, na.rm = TRUE),
    sd = sd(pct_esp_vert_2011ct, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_2011ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)    
## gentrified_2016_2011   1  23284   23284   102.9 <2e-16 ***
## Residuals            688 155714     226                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_2011ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.998 -10.293  -0.814   8.654  50.792 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               41.1580     0.6852   60.06   <2e-16 ***
## gentrified_2016_2011TRUE -12.6588     1.2481  -10.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.04 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1301, Adjusted R-squared:  0.1288 
## F-statistic: 102.9 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.3.2.2 Buffer 250m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_2011b250, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_2011b250, na.rm = TRUE),
    sd = sd(pct_esp_vert_2011b250, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_2011b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)    
## gentrified_2016_2011   1  14776   14776   94.56 <2e-16 ***
## Residuals            688 107511     156                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_2011b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.278  -8.132  -1.213   7.587  40.172 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               38.6277     0.5694  67.841   <2e-16 ***
## gentrified_2016_2011TRUE -10.0844     1.0371  -9.724   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.5 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1208, Adjusted R-squared:  0.1196 
## F-statistic: 94.56 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.3.2.3 Buffer 500m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_2011b500, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_2011b500, na.rm = TRUE),
    sd = sd(pct_esp_vert_2011b500, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_2011b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)    
## gentrified_2016_2011   1  12682   12682   90.12 <2e-16 ***
## Residuals            688  96815     141                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_2011b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.213  -7.841  -1.241   7.464  38.087 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               38.5032     0.5403  71.260   <2e-16 ***
## gentrified_2016_2011TRUE  -9.3426     0.9841  -9.493   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.86 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1158, Adjusted R-squared:  0.1145 
## F-statistic: 90.12 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.3.2.4 Buffer 750m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_2011b750, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_2011b750, na.rm = TRUE),
    sd = sd(pct_esp_vert_2011b750, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_2011b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)    
## gentrified_2016_2011   1  11029   11029   82.95 <2e-16 ***
## Residuals            688  91471     133                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_2011b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_2011b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.141  -7.901  -1.059   7.329  37.659 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               38.4109     0.5252  73.136   <2e-16 ***
## gentrified_2016_2011TRUE  -8.7122     0.9566  -9.108   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.53 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1076, Adjusted R-squared:  0.1063 
## F-statistic: 82.95 on 1 and 688 DF,  p-value: < 2.2e-16

4.1.3.3 Canopy (trees)

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %)

4.1.3.3.1 Census tract level

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_high_2011ct, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_high_2011ct, na.rm = TRUE),
    sd = sd(pct_esp_vert_high_2011ct, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_high_2011ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1   1856  1855.8   20.54 6.89e-06 ***
## Residuals            688  62160    90.3                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_high_2011ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011ct ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.448  -6.230  -1.735   4.263  49.252 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               19.3483     0.4329  44.690  < 2e-16 ***
## gentrified_2016_2011TRUE  -3.5738     0.7886  -4.532 6.89e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.505 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.02899,    Adjusted R-squared:  0.02758 
## F-statistic: 20.54 on 1 and 688 DF,  p-value: 6.886e-06

4.1.3.3.2 Buffer 250m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_high_2011b250, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_high_2011b250, na.rm = TRUE),
    sd = sd(pct_esp_vert_high_2011b250, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_high_2011b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1    786   786.4   13.43 0.000266 ***
## Residuals            688  40275    58.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_high_2011b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b250 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.707  -5.225  -1.519   3.693  31.513 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               17.9168     0.3485  51.411  < 2e-16 ***
## gentrified_2016_2011TRUE  -2.3264     0.6347  -3.665 0.000266 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.651 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01915,    Adjusted R-squared:  0.01773 
## F-statistic: 13.43 on 1 and 688 DF,  p-value: 0.0002662

4.1.3.3.3 Buffer 500m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_high_2011b500, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_high_2011b500, na.rm = TRUE),
    sd = sd(pct_esp_vert_high_2011b500, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_high_2011b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1    632   632.0   12.27 0.000491 ***
## Residuals            688  35449    51.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_high_2011b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b500 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.868  -4.719  -1.365   3.895  31.012 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               17.9576     0.3270  54.924  < 2e-16 ***
## gentrified_2016_2011TRUE  -2.0855     0.5955  -3.502 0.000491 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.178 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01752,    Adjusted R-squared:  0.01609 
## F-statistic: 12.27 on 1 and 688 DF,  p-value: 0.0004913

4.1.3.3.4 Buffer 750m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_high_2011b750, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_high_2011b750, na.rm = TRUE),
    sd = sd(pct_esp_vert_high_2011b750, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_high_2011b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value  Pr(>F)   
## gentrified_2016_2011   1    476   476.4    10.2 0.00147 **
## Residuals            688  32147    46.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_high_2011b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_high_2011b750 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.393  -4.715  -1.352   3.662  28.307 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               17.9127     0.3114  57.532  < 2e-16 ***
## gentrified_2016_2011TRUE  -1.8107     0.5671  -3.193  0.00147 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.836 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0146, Adjusted R-squared:  0.01317 
## F-statistic: 10.19 on 1 and 688 DF,  p-value: 0.001473

5 Association between SES and Urban Interventions at baseline (2011)

Looking at objective #1 | do urban interventions tend to be located in low SES neighborhoods?. We look at \[Urban Intervention_{2011 \to 2016} = f(SES_{2011})\] as well as \[Urban Intervention_{2011 \to 2016} = f(Gentrification_{2011 \to 2016})\]

Here \(Urban Intervention\) means the changes in the urban environment features, such as variation of bike lane length, greenness coverage, etc.

5.1 INTERACT study area

# keep only interact CT
bei_df_aoi <- filter(bei_df, interact_aoi)

5.1.1 UI vs Pampalon | material

5.1.1.1 Bike lane length change

Bike lane ratio to streets (in %)

5.1.1.1.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane_diff.by.street.2011.2016ct, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane_diff.by.street.2011.2016ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016ct ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.442  -3.695  -3.621   0.992  38.452 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.6795     0.2761  13.328   <2e-16 ***
## wSCOREMAT.2011  -1.2310     6.7303  -0.183    0.855    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.194 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  4.862e-05,  Adjusted R-squared:  -0.001405 
## F-statistic: 0.03345 on 1 and 688 DF,  p-value: 0.8549

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016ct ~ wSCOREMAT.2011 + Bike_lane.by.street.2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016ct ~ wSCOREMAT.2011 + 
##     Bike_lane.by.street.2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.601  -4.253  -3.044   1.520  37.493 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 4.62603    0.37517  12.331  < 2e-16 ***
## wSCOREMAT.2011             -4.12185    6.71541  -0.614 0.539558    
## Bike_lane.by.street.2011ct -0.10651    0.02889  -3.687 0.000245 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.129 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01945,    Adjusted R-squared:  0.0166 
## F-statistic: 6.815 on 2 and 687 DF,  p-value: 0.001173

5.1.1.1.2 Buffers 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane_diff.by.street.2011.2016b250, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane_diff.by.street.2011.2016b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b250 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.425  -3.455  -2.437   1.626  27.938 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.4854     0.2095  16.637   <2e-16 ***
## wSCOREMAT.2011  -5.3511     5.1073  -1.048    0.295    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.459 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.001593,   Adjusted R-squared:  0.0001418 
## F-statistic: 1.098 on 1 and 688 DF,  p-value: 0.2951

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b250 ~ wSCOREMAT.2011 + Bike_lane.by.street.2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b250 ~ wSCOREMAT.2011 + 
##     Bike_lane.by.street.2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.103  -3.597  -1.898   1.815  26.261 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    5.20649    0.33995  15.315  < 2e-16 ***
## wSCOREMAT.2011               -11.81363    5.07233  -2.329   0.0201 *  
## Bike_lane.by.street.2011b250  -0.19320    0.03054  -6.325 4.55e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.311 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05654,    Adjusted R-squared:  0.05379 
## F-statistic: 20.59 on 2 and 687 DF,  p-value: 2.078e-09

5.1.1.1.3 Buffers 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane_diff.by.street.2011.2016b500, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane_diff.by.street.2011.2016b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b500 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.367 -3.310 -1.827  2.005 18.540 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.4779     0.1803   19.29   <2e-16 ***
## wSCOREMAT.2011  -6.4606     4.3960   -1.47    0.142    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.699 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.00313,    Adjusted R-squared:  0.001681 
## F-statistic:  2.16 on 1 and 688 DF,  p-value: 0.1421

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b500 ~ wSCOREMAT.2011 + Bike_lane.by.street.2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b500 ~ wSCOREMAT.2011 + 
##     Bike_lane.by.street.2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.405 -3.254 -1.529  2.070 18.630 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    5.15019    0.32914  15.648  < 2e-16 ***
## wSCOREMAT.2011               -12.25032    4.39472  -2.788  0.00546 ** 
## Bike_lane.by.street.2011b500  -0.18645    0.03102  -6.011 2.99e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.583 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05294,    Adjusted R-squared:  0.05018 
## F-statistic:  19.2 on 2 and 687 DF,  p-value: 7.689e-09

5.1.1.1.4 Buffers 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane_diff.by.street.2011.2016b750, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane_diff.by.street.2011.2016b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b750 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.565 -3.039 -1.415  1.939 18.956 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.4705     0.1622   21.40   <2e-16 ***
## wSCOREMAT.2011  -5.5739     3.9538   -1.41    0.159    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.226 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.00288,    Adjusted R-squared:  0.001431 
## F-statistic: 1.987 on 1 and 688 DF,  p-value: 0.1591

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b750 ~ wSCOREMAT.2011 + Bike_lane.by.street.2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b750 ~ wSCOREMAT.2011 + 
##     Bike_lane.by.street.2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.941 -3.024 -1.105  1.910 18.018 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    5.12804    0.32188  15.931  < 2e-16 ***
## wSCOREMAT.2011               -10.52276    3.94933  -2.664  0.00789 ** 
## Bike_lane.by.street.2011b750  -0.18420    0.03114  -5.914 5.25e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.126 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05119,    Adjusted R-squared:  0.04843 
## F-statistic: 18.53 on 2 and 687 DF,  p-value: 1.449e-08

5.1.1.2 Canopy change

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %)

5.1.1.2.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_2011.2017ct, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_2011.2017ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017ct ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.165  -1.944  -0.273   1.524  37.029 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      4.1619     0.1313  31.687   <2e-16 ***
## wSCOREMAT.2011   1.7219     3.2020   0.538    0.591    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.423 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0004201,  Adjusted R-squared:  -0.001033 
## F-statistic: 0.2892 on 1 and 688 DF,  p-value: 0.5909

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017ct ~ wSCOREMAT.2011 + pct_esp_vert_2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017ct ~ wSCOREMAT.2011 + 
##     pct_esp_vert_2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.137  -1.684  -0.296   1.369  36.968 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          5.936893   0.338723  17.527  < 2e-16 ***
## wSCOREMAT.2011      -3.393906   3.259702  -1.041    0.298    
## pct_esp_vert_2011ct -0.046757   0.008256  -5.663 2.18e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.348 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04501,    Adjusted R-squared:  0.04223 
## F-statistic: 16.19 on 2 and 687 DF,  p-value: 1.35e-07

5.1.1.2.2 Buffer 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_2011.2017b250, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_2011.2017b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b250 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.827 -1.729 -0.027  1.341 34.699 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.9602     0.1125  35.192   <2e-16 ***
## wSCOREMAT.2011  -1.7145     2.7434  -0.625    0.532    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.932 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0005674,  Adjusted R-squared:  -0.0008853 
## F-statistic: 0.3906 on 1 and 688 DF,  p-value: 0.5322

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b250 ~ wSCOREMAT.2011 + pct_esp_vert_2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b250 ~ wSCOREMAT.2011 + 
##     pct_esp_vert_2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.958 -1.514 -0.185  1.094 34.802 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            5.95375    0.32597  18.265  < 2e-16 ***
## wSCOREMAT.2011        -6.73103    2.77465  -2.426   0.0155 *  
## pct_esp_vert_2011b250 -0.05524    0.00851  -6.492 1.63e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.849 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05833,    Adjusted R-squared:  0.05559 
## F-statistic: 21.28 on 2 and 687 DF,  p-value: 1.082e-09

5.1.1.2.3 Buffer 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_2011.2017b500, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_2011.2017b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b500 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.614 -1.542 -0.055  1.293 32.079 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.8940     0.1049  37.116   <2e-16 ***
## wSCOREMAT.2011  -1.2350     2.5576  -0.483    0.629    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.734 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0003388,  Adjusted R-squared:  -0.001114 
## F-statistic: 0.2332 on 1 and 688 DF,  p-value: 0.6293

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b500 ~ wSCOREMAT.2011 + pct_esp_vert_2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b500 ~ wSCOREMAT.2011 + 
##     pct_esp_vert_2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.817 -1.482 -0.186  1.122 32.262 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            5.75759    0.32003  17.991  < 2e-16 ***
## wSCOREMAT.2011        -5.52570    2.58791  -2.135   0.0331 *  
## pct_esp_vert_2011b500 -0.05156    0.00839  -6.145 1.35e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.664 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05242,    Adjusted R-squared:  0.04967 
## F-statistic:    19 on 2 and 687 DF,  p-value: 9.264e-09

5.1.1.2.4 Buffer 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_2011.2017b750, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_2011.2017b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b750 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6691 -1.5113 -0.0059  1.2007 29.1958 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.8425     0.1008  38.102   <2e-16 ***
## wSCOREMAT.2011  -1.7355     2.4585  -0.706     0.48    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.628 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0007238,  Adjusted R-squared:  -0.0007286 
## F-statistic: 0.4983 on 1 and 688 DF,  p-value: 0.4805

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b750 ~ wSCOREMAT.2011 + pct_esp_vert_2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b750 ~ wSCOREMAT.2011 + 
##     pct_esp_vert_2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4480 -1.4348 -0.1886  1.1088 29.4180 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            5.619194   0.317052  17.723  < 2e-16 ***
## wSCOREMAT.2011        -5.470665   2.482528  -2.204   0.0279 *  
## pct_esp_vert_2011b750 -0.049066   0.008323  -5.895 5.86e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.566 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04884,    Adjusted R-squared:  0.04607 
## F-statistic: 17.64 on 2 and 687 DF,  p-value: 3.385e-08

5.1.1.3 Canopy (trees) change

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %)

5.1.1.3.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_high_2011.2017ct, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_high_2011.2017ct ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017ct ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8508 -1.4234 -0.3206  1.1810 11.3606 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.4314     0.0846  40.560  < 2e-16 ***
## wSCOREMAT.2011 -11.9113     2.0625  -5.775 1.16e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.205 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04624,    Adjusted R-squared:  0.04485 
## F-statistic: 33.35 on 1 and 688 DF,  p-value: 1.164e-08

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017ct ~ wSCOREMAT.2011 + pct_esp_vert_high_2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017ct ~ wSCOREMAT.2011 + 
##     pct_esp_vert_high_2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.857 -1.451 -0.343  1.155 10.810 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               2.945265   0.197553  14.909  < 2e-16 ***
## wSCOREMAT.2011           -9.410419   2.249397  -4.184 3.24e-05 ***
## pct_esp_vert_high_2011ct  0.025850   0.009502   2.720  0.00668 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.194 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0564, Adjusted R-squared:  0.05365 
## F-statistic: 20.53 on 2 and 687 DF,  p-value: 2.184e-09

5.1.1.3.2 Buffer 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_high_2011.2017b250, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_high_2011.2017b250 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b250 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9628 -1.1429 -0.2624  0.8870  6.8083 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.15716    0.06809  46.369  < 2e-16 ***
## wSCOREMAT.2011 -12.18518    1.65991  -7.341 6.01e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.774 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.07264,    Adjusted R-squared:  0.07129 
## F-statistic: 53.89 on 1 and 688 DF,  p-value: 6.007e-13

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b250 ~ wSCOREMAT.2011 + pct_esp_vert_high_2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b250 ~ wSCOREMAT.2011 + 
##     pct_esp_vert_high_2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4856 -1.1112 -0.2855  0.8853  6.7932 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  2.694806   0.183921  14.652  < 2e-16 ***
## wSCOREMAT.2011             -10.043738   1.832306  -5.481 5.93e-08 ***
## pct_esp_vert_high_2011b250   0.026190   0.009685   2.704  0.00702 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.766 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0824, Adjusted R-squared:  0.07973 
## F-statistic: 30.85 on 2 and 687 DF,  p-value: 1.482e-13

5.1.1.3.3 Buffer 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_high_2011.2017b500, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_high_2011.2017b500 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b500 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0396 -1.0949 -0.2033  0.7381  7.3251 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.07086    0.06356  48.316  < 2e-16 ***
## wSCOREMAT.2011 -10.94933    1.54945  -7.067  3.9e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.656 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06767,    Adjusted R-squared:  0.06632 
## F-statistic: 49.94 on 1 and 688 DF,  p-value: 3.902e-12

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b500 ~ wSCOREMAT.2011 + pct_esp_vert_high_2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b500 ~ wSCOREMAT.2011 + 
##     pct_esp_vert_high_2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8059 -1.0936 -0.1912  0.6904  7.2757 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 2.558316   0.182249  14.037  < 2e-16 ***
## wSCOREMAT.2011             -8.727890   1.709430  -5.106 4.27e-07 ***
## pct_esp_vert_high_2011b500  0.028895   0.009637   2.998  0.00281 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.647 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.07971,    Adjusted R-squared:  0.07703 
## F-statistic: 29.75 on 2 and 687 DF,  p-value: 4.051e-13

5.1.1.3.4 Buffer 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_high_2011.2017b750, x=wSCOREMAT.2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_high_2011.2017b750 ~ wSCOREMAT.2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b750 ~ wSCOREMAT.2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9569 -0.9961 -0.2127  0.7393  6.2521 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.99802    0.05853  51.218  < 2e-16 ***
## wSCOREMAT.2011 -10.33948    1.42699  -7.246 1.16e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0709, Adjusted R-squared:  0.06955 
## F-statistic:  52.5 on 1 and 688 DF,  p-value: 1.158e-12

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b750 ~ wSCOREMAT.2011 + pct_esp_vert_high_2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b750 ~ wSCOREMAT.2011 + 
##     pct_esp_vert_high_2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7591 -0.9842 -0.2096  0.6951  6.2118 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 2.439850   0.175338  13.915  < 2e-16 ***
## wSCOREMAT.2011             -8.045919   1.571026  -5.121 3.94e-07 ***
## pct_esp_vert_high_2011b750  0.031433   0.009316   3.374 0.000783 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.514 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.08604,    Adjusted R-squared:  0.08338 
## F-statistic: 32.34 on 2 and 687 DF,  p-value: 3.786e-14

5.1.2 UI vs Visible minority

5.1.2.1 Bike lane length change

Bike lane ratio to streets (in %). With or without controlling for UC in 2011.

5.1.2.1.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane_diff.by.street.2011.2016ct, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane_diff.by.street.2011.2016ct ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016ct ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.758 -3.975 -3.059  0.859 38.865 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4.85228    0.51159   9.485  < 2e-16 ***
## vis_minority_2011 -0.04599    0.01689  -2.723  0.00663 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.156 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01066,    Adjusted R-squared:  0.009226 
## F-statistic: 7.416 on 1 and 688 DF,  p-value: 0.00663

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016ct ~ vis_minority_2011 + Bike_lane.by.street.2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016ct ~ vis_minority_2011 + 
##     Bike_lane.by.street.2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.788 -4.145 -2.856  1.391 37.922 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 6.17860    0.59938  10.308  < 2e-16 ***
## vis_minority_2011          -0.05708    0.01691  -3.375 0.000779 ***
## Bike_lane.by.street.2011ct -0.11971    0.02903  -4.123  4.2e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.074 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.03455,    Adjusted R-squared:  0.03174 
## F-statistic: 12.29 on 2 and 687 DF,  p-value: 5.679e-06

5.1.2.1.2 Buffers 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane_diff.by.street.2011.2016b250, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane_diff.by.street.2011.2016b250 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b250 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.625  -3.447  -2.147   1.972  27.232 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4.86710    0.38548  12.626  < 2e-16 ***
## vis_minority_2011 -0.05521    0.01273  -4.338 1.65e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.392 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.02662,    Adjusted R-squared:  0.02521 
## F-statistic: 18.82 on 1 and 688 DF,  p-value: 1.655e-05

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b250 ~ vis_minority_2011 + Bike_lane.by.street.2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b250 ~ vis_minority_2011 + 
##     Bike_lane.by.street.2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.902  -3.524  -1.576   2.039  24.834 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   7.53777    0.51119  14.745  < 2e-16 ***
## vis_minority_2011            -0.08075    0.01269  -6.364 3.60e-10 ***
## Bike_lane.by.street.2011b250 -0.23135    0.03051  -7.584 1.09e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.183 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1018, Adjusted R-squared:  0.0992 
## F-statistic: 38.94 on 2 and 687 DF,  p-value: < 2.2e-16

5.1.2.1.3 Buffers 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane_diff.by.street.2011.2016b500, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane_diff.by.street.2011.2016b500 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b500 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.704 -3.241 -1.655  2.003 18.311 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4.59890    0.33235  13.837  < 2e-16 ***
## vis_minority_2011 -0.04575    0.01097  -4.169 3.45e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.649 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.02464,    Adjusted R-squared:  0.02323 
## F-statistic: 17.38 on 1 and 688 DF,  p-value: 3.446e-05

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b500 ~ vis_minority_2011 + Bike_lane.by.street.2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b500 ~ vis_minority_2011 + 
##     Bike_lane.by.street.2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.221 -3.234 -1.230  1.929 18.372 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   7.19708    0.48166  14.942  < 2e-16 ***
## vis_minority_2011            -0.06964    0.01109  -6.280 6.01e-10 ***
## Bike_lane.by.street.2011b500 -0.22591    0.03125  -7.229 1.30e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.485 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.09358,    Adjusted R-squared:  0.09094 
## F-statistic: 35.46 on 2 and 687 DF,  p-value: 2.199e-15

5.1.2.1.4 Buffers 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=Bike_lane_diff.by.street.2011.2016b750, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(Bike_lane_diff.by.street.2011.2016b750 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b750 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.988 -3.000 -1.417  1.921 18.327 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4.45024    0.29897  14.885  < 2e-16 ***
## vis_minority_2011 -0.04016    0.00987  -4.068 5.29e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.182 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.02349,    Adjusted R-squared:  0.02207 
## F-statistic: 16.55 on 1 and 688 DF,  p-value: 5.286e-05

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b750 ~ vis_minority_2011 + Bike_lane.by.street.2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b750 ~ vis_minority_2011 + 
##     Bike_lane.by.street.2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.973 -2.915 -1.020  1.842 16.865 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   7.07360    0.46138  15.331  < 2e-16 ***
## vis_minority_2011            -0.06332    0.01003  -6.310 5.01e-10 ***
## Bike_lane.by.street.2011b750 -0.22946    0.03151  -7.282 9.02e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.032 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.09347,    Adjusted R-squared:  0.09083 
## F-statistic: 35.42 on 2 and 687 DF,  p-value: 2.296e-15

5.1.2.2 Canopy change

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %)

5.1.2.2.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_2011.2017ct, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_2011.2017ct ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017ct ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.093  -1.938  -0.199   1.553  36.892 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4.257351   0.244837   17.39   <2e-16 ***
## vis_minority_2011 -0.004040   0.008083   -0.50    0.617    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.425 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.000363,   Adjusted R-squared:  -0.00109 
## F-statistic: 0.2498 on 1 and 688 DF,  p-value: 0.6174

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017ct ~ vis_minority_2011 + pct_esp_vert_2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017ct ~ vis_minority_2011 + 
##     pct_esp_vert_2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.326  -1.777  -0.310   1.361  36.775 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          6.230912   0.423236  14.722  < 2e-16 ***
## vis_minority_2011   -0.014146   0.008107  -1.745   0.0814 .  
## pct_esp_vert_2011ct -0.045913   0.008118  -5.656 2.28e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.35 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04484,    Adjusted R-squared:  0.04205 
## F-statistic: 16.12 on 2 and 687 DF,  p-value: 1.435e-07

5.1.2.2.2 Buffer 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_2011.2017b250, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_2011.2017b250 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b250 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.745 -1.678 -0.043  1.310 34.485 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4.253273   0.209119  20.339   <2e-16 ***
## vis_minority_2011 -0.011863   0.006904  -1.718   0.0862 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.925 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.004273,   Adjusted R-squared:  0.002825 
## F-statistic: 2.952 on 1 and 688 DF,  p-value: 0.08621

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b250 ~ vis_minority_2011 + pct_esp_vert_2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b250 ~ vis_minority_2011 + 
##     pct_esp_vert_2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.593 -1.584 -0.145  1.141 34.553 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            6.400991   0.387333  16.526  < 2e-16 ***
## vis_minority_2011     -0.020710   0.006842  -3.027  0.00256 ** 
## pct_esp_vert_2011b250 -0.053976   0.008289  -6.512 1.43e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.841 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06216,    Adjusted R-squared:  0.05943 
## F-statistic: 22.77 on 2 and 687 DF,  p-value: 2.668e-10

5.1.2.2.3 Buffer 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_2011.2017b500, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_2011.2017b500 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b500 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.352 -1.594 -0.054  1.310 31.816 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4.225329   0.194825  21.688   <2e-16 ***
## vis_minority_2011 -0.013070   0.006432  -2.032   0.0425 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.725 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.005966,   Adjusted R-squared:  0.004521 
## F-statistic: 4.129 on 1 and 688 DF,  p-value: 0.04253

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b500 ~ vis_minority_2011 + pct_esp_vert_2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b500 ~ vis_minority_2011 + 
##     pct_esp_vert_2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.510 -1.493 -0.199  1.206 31.983 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            6.241340   0.371498  16.800  < 2e-16 ***
## vis_minority_2011     -0.020193   0.006359  -3.176  0.00156 ** 
## pct_esp_vert_2011b500 -0.051374   0.008142  -6.310    5e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.651 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06042,    Adjusted R-squared:  0.05768 
## F-statistic: 22.09 on 2 and 687 DF,  p-value: 5.046e-10

5.1.2.2.4 Buffer 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_2011.2017b750, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_2011.2017b750 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b750 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6430 -1.5699 -0.0092  1.2938 28.8848 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4.242668   0.186996   22.69   <2e-16 ***
## vis_minority_2011 -0.015868   0.006174   -2.57   0.0104 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.616 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.009511,   Adjusted R-squared:  0.008071 
## F-statistic: 6.606 on 1 and 688 DF,  p-value: 0.01037

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b750 ~ vis_minority_2011 + pct_esp_vert_2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b750 ~ vis_minority_2011 + 
##     pct_esp_vert_2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6188 -1.4098 -0.1884  1.1195 29.1015 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            6.148210   0.362749  16.949  < 2e-16 ***
## vis_minority_2011     -0.021772   0.006096  -3.571  0.00038 ***
## pct_esp_vert_2011b750 -0.049020   0.008068  -6.076 2.04e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.55 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06002,    Adjusted R-squared:  0.05729 
## F-statistic: 21.93 on 2 and 687 DF,  p-value: 5.83e-10

5.1.2.3 Canopy (trees) change

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %)

5.1.2.3.1 Census tract level

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_high_2011.2017ct, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_high_2011.2017ct ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017ct ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7482 -1.5184 -0.2624  1.1531 10.7915 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.990009   0.159253  25.055  < 2e-16 ***
## vis_minority_2011 -0.024453   0.005258  -4.651 3.96e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.228 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.03048,    Adjusted R-squared:  0.02907 
## F-statistic: 21.63 on 1 and 688 DF,  p-value: 3.964e-06

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017ct ~ vis_minority_2011 + pct_esp_vert_high_2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017ct ~ vis_minority_2011 + 
##     pct_esp_vert_high_2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7699 -1.4988 -0.3498  1.1095 10.2885 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               3.200485   0.251171  12.742  < 2e-16 ***
## vis_minority_2011        -0.019413   0.005348  -3.630 0.000305 ***
## pct_esp_vert_high_2011ct  0.036140   0.008956   4.035 6.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.203 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05293,    Adjusted R-squared:  0.05017 
## F-statistic:  19.2 on 2 and 687 DF,  p-value: 7.708e-09

5.1.2.3.2 Buffer 250m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_high_2011.2017b250, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_high_2011.2017b250 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b250 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0818 -1.2094 -0.2710  0.9444  6.9013 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.657494   0.129117  28.327  < 2e-16 ***
## vis_minority_2011 -0.021890   0.004263  -5.135 3.67e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.806 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.03691,    Adjusted R-squared:  0.03551 
## F-statistic: 26.37 on 1 and 688 DF,  p-value: 3.673e-07

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b250 ~ vis_minority_2011 + pct_esp_vert_high_2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b250 ~ vis_minority_2011 + 
##     pct_esp_vert_high_2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1679 -1.1641 -0.2817  0.8024  6.8128 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 2.843690   0.217583  13.069  < 2e-16 ***
## vis_minority_2011          -0.017914   0.004289  -4.177 3.34e-05 ***
## pct_esp_vert_high_2011b250  0.041350   0.008967   4.611 4.78e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.78 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06583,    Adjusted R-squared:  0.06311 
## F-statistic:  24.2 on 2 and 687 DF,  p-value: 6.948e-11

5.1.2.3.3 Buffer 500m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_high_2011.2017b500, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_high_2011.2017b500 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b500 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2936 -1.1328 -0.2707  0.8096  7.3055 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.474949   0.120735  28.782  < 2e-16 ***
## vis_minority_2011 -0.017844   0.003986  -4.477 8.87e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.689 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0283, Adjusted R-squared:  0.02689 
## F-statistic: 20.04 on 1 and 688 DF,  p-value: 8.873e-06

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b500 ~ vis_minority_2011 + pct_esp_vert_high_2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b500 ~ vis_minority_2011 + 
##     pct_esp_vert_high_2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4292 -1.0788 -0.2186  0.7476  7.1924 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 2.624779   0.207495  12.650  < 2e-16 ***
## vis_minority_2011          -0.014596   0.003972  -3.675 0.000257 ***
## pct_esp_vert_high_2011b500  0.044255   0.008859   4.995 7.46e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.66 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06236,    Adjusted R-squared:  0.05963 
## F-statistic: 22.84 on 2 and 687 DF,  p-value: 2.481e-10

5.1.2.3.4 Buffer 750m

ggplot(units::drop_units(bei_df_aoi), aes(y=pct_esp_vert_diff_high_2011.2017b750, x=vis_minority_2011)) +
  geom_point() +
  geom_smooth(method=lm)

res.lm <- lm(pct_esp_vert_diff_high_2011.2017b750 ~ vis_minority_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b750 ~ vis_minority_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1753 -1.0349 -0.2860  0.7644  6.2564 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.347398   0.111550  30.008  < 2e-16 ***
## vis_minority_2011 -0.015569   0.003683  -4.228 2.68e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.56 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.02532,    Adjusted R-squared:  0.0239 
## F-statistic: 17.87 on 1 and 688 DF,  p-value: 2.681e-05

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b750 ~ vis_minority_2011 + pct_esp_vert_high_2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b750 ~ vis_minority_2011 + 
##     pct_esp_vert_high_2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3455 -1.0266 -0.2463  0.6975  6.1498 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 2.458366   0.194986  12.608  < 2e-16 ***
## vis_minority_2011          -0.012756   0.003643  -3.502 0.000493 ***
## pct_esp_vert_high_2011b750  0.047037   0.008545   5.505 5.23e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.528 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06649,    Adjusted R-squared:  0.06378 
## F-statistic: 24.47 on 2 and 687 DF,  p-value: 5.437e-11

5.1.3 UI vs gentrified CT

Gentrified CT between 2011 and 2016

5.1.3.1 Bike lane length change

Bike lane ratio to streets (in %)

5.1.3.1.1 Census tract level

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=Bike_lane_diff.by.street.2011.2016ct, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(Bike_lane_diff.by.street.2011.2016ct, na.rm = TRUE),
    sd = sd(Bike_lane_diff.by.street.2011.2016ct, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(Bike_lane_diff.by.street.2011.2016ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1    843   843.4   16.69 4.92e-05 ***
## Residuals            688  34767    50.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(Bike_lane_diff.by.street.2011.2016ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016ct ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.724 -2.947 -2.947  0.662 39.191 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                2.9468     0.3238   9.101  < 2e-16 ***
## gentrified_2016_2011TRUE   2.4092     0.5897   4.085 4.92e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.109 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.02368,    Adjusted R-squared:  0.02226 
## F-statistic: 16.69 on 1 and 688 DF,  p-value: 4.921e-05

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016ct ~ gentrified_2016_2011 + Bike_lane.by.street.2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016ct ~ gentrified_2016_2011 + 
##     Bike_lane.by.street.2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.898 -3.850 -2.671  1.060 38.291 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 3.84970    0.40696   9.460  < 2e-16 ***
## gentrified_2016_2011TRUE    2.39022    0.58467   4.088 4.86e-05 ***
## Bike_lane.by.street.2011ct -0.10307    0.02856  -3.609 0.000329 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.047 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04185,    Adjusted R-squared:  0.03906 
## F-statistic:    15 on 2 and 687 DF,  p-value: 4.19e-07

5.1.3.1.2 Buffers 250m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=Bike_lane_diff.by.street.2011.2016b250, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(Bike_lane_diff.by.street.2011.2016b250, na.rm = TRUE),
    sd = sd(Bike_lane_diff.by.street.2011.2016b250, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(Bike_lane_diff.by.street.2011.2016b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1   1124  1123.9    39.8 5.03e-10 ***
## Residuals            688  19426    28.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b250 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.617  -2.613  -2.285   1.770  26.114 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                2.6134     0.2420  10.798  < 2e-16 ***
## gentrified_2016_2011TRUE   2.7811     0.4408   6.309 5.03e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.314 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05469,    Adjusted R-squared:  0.05332 
## F-statistic:  39.8 on 1 and 688 DF,  p-value: 5.027e-10

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b250 ~ gentrified_2016_2011 + Bike_lane.by.street.2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b250 ~ gentrified_2016_2011 + 
##     Bike_lane.by.street.2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.345  -3.320  -1.535   1.761  24.816 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   4.18653    0.34725  12.056  < 2e-16 ***
## gentrified_2016_2011TRUE      2.80059    0.42942   6.522 1.35e-10 ***
## Bike_lane.by.street.2011b250 -0.18123    0.02937  -6.170 1.16e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.176 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1043, Adjusted R-squared:  0.1017 
## F-statistic: 40.01 on 2 and 687 DF,  p-value: < 2.2e-16

5.1.3.1.3 Buffers 500m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=Bike_lane_diff.by.street.2011.2016b500, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(Bike_lane_diff.by.street.2011.2016b500, na.rm = TRUE),
    sd = sd(Bike_lane_diff.by.street.2011.2016b500, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(Bike_lane_diff.by.street.2011.2016b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1    965   965.1    46.5 2.01e-11 ***
## Residuals            688  14280    20.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b500 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.615 -2.649 -1.577  1.713 18.848 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                2.6491     0.2075  12.766  < 2e-16 ***
## gentrified_2016_2011TRUE   2.5772     0.3779   6.819 2.01e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.556 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06331,    Adjusted R-squared:  0.06195 
## F-statistic:  46.5 on 1 and 688 DF,  p-value: 2.008e-11

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b500 ~ gentrified_2016_2011 + Bike_lane.by.street.2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b500 ~ gentrified_2016_2011 + 
##     Bike_lane.by.street.2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.993 -2.865 -1.310  1.756 18.065 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   4.15807    0.32708  12.713  < 2e-16 ***
## gentrified_2016_2011TRUE      2.64460    0.36924   7.162 2.05e-12 ***
## Bike_lane.by.street.2011b500 -0.17399    0.02961  -5.877 6.52e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.449 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1081, Adjusted R-squared:  0.1055 
## F-statistic: 41.65 on 2 and 687 DF,  p-value: < 2.2e-16

5.1.3.1.4 Buffers 750m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=Bike_lane_diff.by.street.2011.2016b750, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(Bike_lane_diff.by.street.2011.2016b750, na.rm = TRUE),
    sd = sd(Bike_lane_diff.by.street.2011.2016b750, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(Bike_lane_diff.by.street.2011.2016b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1    759   758.6   45.14 3.85e-11 ***
## Residuals            688  11563    16.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b750 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.893 -2.732 -1.275  1.702 17.528 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                2.7319     0.1867  14.630  < 2e-16 ***
## gentrified_2016_2011TRUE   2.2850     0.3401   6.719 3.85e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.1 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.06157,    Adjusted R-squared:  0.06021 
## F-statistic: 45.14 on 1 and 688 DF,  p-value: 3.848e-11

# Accounting for UC in 2011
res.lm <- lm(Bike_lane_diff.by.street.2011.2016b750 ~ gentrified_2016_2011 + Bike_lane.by.street.2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = Bike_lane_diff.by.street.2011.2016b750 ~ gentrified_2016_2011 + 
##     Bike_lane.by.street.2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1464 -2.7321 -0.9453  1.7715 16.6912 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   4.26493    0.31573  13.508  < 2e-16 ***
## gentrified_2016_2011TRUE      2.37858    0.33229   7.158 2.10e-12 ***
## Bike_lane.by.street.2011b750 -0.17651    0.02969  -5.946 4.38e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.001 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1075, Adjusted R-squared:  0.1049 
## F-statistic: 41.37 on 2 and 687 DF,  p-value: < 2.2e-16

5.1.3.2 Canopy change

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %)

5.1.3.2.1 Census tract level

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_diff_2011.2017ct, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_diff_2011.2017ct, na.rm = TRUE),
    sd = sd(pct_esp_vert_diff_2011.2017ct, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_diff_2011.2017ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value  Pr(>F)   
## gentrified_2016_2011   1     97   96.85   8.355 0.00397 **
## Residuals            688   7975   11.59                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_diff_2011.2017ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017ct ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.644  -1.868  -0.258   1.425  37.232 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                3.9077     0.1551   25.20  < 2e-16 ***
## gentrified_2016_2011TRUE   0.8164     0.2825    2.89  0.00397 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.405 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.012,  Adjusted R-squared:  0.01056 
## F-statistic: 8.355 on 1 and 688 DF,  p-value: 0.003968

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017ct ~ gentrified_2016_2011 + pct_esp_vert_2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017ct ~ gentrified_2016_2011 + 
##     pct_esp_vert_2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.455  -1.704  -0.312   1.382  37.183 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.535249   0.381819  14.497  < 2e-16 ***
## gentrified_2016_2011TRUE  0.315837   0.298396   1.058     0.29    
## pct_esp_vert_2011ct      -0.039545   0.008502  -4.651 3.95e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.355 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04216,    Adjusted R-squared:  0.03938 
## F-statistic: 15.12 on 2 and 687 DF,  p-value: 3.746e-07

5.1.3.2.2 Buffer 250m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_diff_2011.2017b250, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_diff_2011.2017b250, na.rm = TRUE),
    sd = sd(pct_esp_vert_diff_2011.2017b250, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_diff_2011.2017b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)  
## gentrified_2016_2011   1     57   56.50   6.639 0.0102 *
## Residuals            688   5855    8.51                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_diff_2011.2017b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b250 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.651 -1.681 -0.051  1.254 34.949 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                3.7612     0.1329  28.305   <2e-16 ***
## gentrified_2016_2011TRUE   0.6236     0.2420   2.577   0.0102 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.917 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.009557,   Adjusted R-squared:  0.008118 
## F-statistic: 6.639 on 1 and 688 DF,  p-value: 0.01018

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b250 ~ gentrified_2016_2011 + pct_esp_vert_2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b250 ~ gentrified_2016_2011 + 
##     pct_esp_vert_2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.849 -1.516 -0.208  1.076 35.045 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.58548    0.36112  15.467  < 2e-16 ***
## gentrified_2016_2011TRUE  0.14733    0.25296   0.582     0.56    
## pct_esp_vert_2011b250    -0.04723    0.00872  -5.416 8.42e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.859 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.05012,    Adjusted R-squared:  0.04735 
## F-statistic: 18.12 on 2 and 687 DF,  p-value: 2.134e-08

5.1.3.2.3 Buffer 500m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_diff_2011.2017b500, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_diff_2011.2017b500, na.rm = TRUE),
    sd = sd(pct_esp_vert_diff_2011.2017b500, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_diff_2011.2017b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value Pr(>F)   
## gentrified_2016_2011   1     59   58.84   7.967 0.0049 **
## Residuals            688   5081    7.39                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_diff_2011.2017b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b500 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.458 -1.598 -0.028  1.220 32.312 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                3.6984     0.1238  29.878   <2e-16 ***
## gentrified_2016_2011TRUE   0.6364     0.2255   2.823   0.0049 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.718 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01145,    Adjusted R-squared:  0.01001 
## F-statistic: 7.967 on 1 and 688 DF,  p-value: 0.004902

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b500 ~ gentrified_2016_2011 + pct_esp_vert_2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b500 ~ gentrified_2016_2011 + 
##     pct_esp_vert_2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.001 -1.447 -0.232  1.035 32.473 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.392180   0.351937  15.321  < 2e-16 ***
## gentrified_2016_2011TRUE  0.225373   0.235477   0.957    0.339    
## pct_esp_vert_2011b500    -0.043991   0.008578  -5.128  3.8e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.669 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0479, Adjusted R-squared:  0.04512 
## F-statistic: 17.28 on 2 and 687 DF,  p-value: 4.764e-08

5.1.3.2.4 Buffer 750m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_diff_2011.2017b750, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_diff_2011.2017b750, na.rm = TRUE),
    sd = sd(pct_esp_vert_diff_2011.2017b750, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_diff_2011.2017b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value  Pr(>F)   
## gentrified_2016_2011   1     53   53.32   7.808 0.00535 **
## Residuals            688   4699    6.83                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_diff_2011.2017b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b750 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8790 -1.5176 -0.0082  1.1328 29.4368 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                3.6532     0.1190  30.690  < 2e-16 ***
## gentrified_2016_2011TRUE   0.6058     0.2168   2.794  0.00535 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.613 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01122,    Adjusted R-squared:  0.009784 
## F-statistic: 7.808 on 1 and 688 DF,  p-value: 0.005348

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_2011.2017b750 ~ gentrified_2016_2011 + pct_esp_vert_2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_2011.2017b750 ~ gentrified_2016_2011 + 
##     pct_esp_vert_2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.548 -1.391 -0.188  1.068 29.627 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.243475   0.346926  15.114  < 2e-16 ***
## gentrified_2016_2011TRUE  0.245108   0.225808   1.085    0.278    
## pct_esp_vert_2011b750    -0.041401   0.008502  -4.870 1.39e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.571 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.04421,    Adjusted R-squared:  0.04143 
## F-statistic: 15.89 on 2 and 687 DF,  p-value: 1.795e-07

5.1.3.3 Canopy (trees) change

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %)

5.1.3.3.1 Census tract level

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_diff_high_2011.2017ct, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_diff_high_2011.2017ct, na.rm = TRUE),
    sd = sd(pct_esp_vert_diff_high_2011.2017ct, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_diff_high_2011.2017ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value  Pr(>F)    
## gentrified_2016_2011   1    106  106.33   21.42 4.4e-06 ***
## Residuals            688   3415    4.96                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_diff_high_2011.2017ct ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017ct ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6852 -1.4752 -0.3852  1.1357 10.1048 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                3.1052     0.1015  30.600  < 2e-16 ***
## gentrified_2016_2011TRUE   0.8555     0.1848   4.628  4.4e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.228 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0302, Adjusted R-squared:  0.02879 
## F-statistic: 21.42 on 1 and 688 DF,  p-value: 4.404e-06

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017ct ~ gentrified_2016_2011 + pct_esp_vert_high_2011ct, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017ct ~ gentrified_2016_2011 + 
##     pct_esp_vert_high_2011ct, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6023 -1.4148 -0.3331  1.0547  9.9179 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              2.095551   0.195587  10.714  < 2e-16 ***
## gentrified_2016_2011TRUE 1.041958   0.182991   5.694 1.84e-08 ***
## pct_esp_vert_high_2011ct 0.052183   0.008718   5.986 3.47e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.174 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.07827,    Adjusted R-squared:  0.07558 
## F-statistic: 29.17 on 2 and 687 DF,  p-value: 6.948e-13

5.1.3.3.2 Buffer 250m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_diff_high_2011.2017b250, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_diff_high_2011.2017b250, na.rm = TRUE),
    sd = sd(pct_esp_vert_diff_high_2011.2017b250, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_diff_high_2011.2017b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1   78.7   78.70   24.05 1.17e-06 ***
## Residuals            688 2251.4    3.27                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b250 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b250 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8444 -1.2019 -0.3344  0.7941  6.9256 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                2.8744     0.0824  34.885  < 2e-16 ***
## gentrified_2016_2011TRUE   0.7360     0.1501   4.904 1.17e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.809 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.03378,    Adjusted R-squared:  0.03237 
## F-statistic: 24.05 on 1 and 688 DF,  p-value: 1.173e-06

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b250 ~ gentrified_2016_2011 + pct_esp_vert_high_2011b250, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b250 ~ gentrified_2016_2011 + 
##     pct_esp_vert_high_2011b250, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0087 -1.1333 -0.3036  0.7533  6.7254 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.870862   0.176274  10.613  < 2e-16 ***
## gentrified_2016_2011TRUE   0.866275   0.147325   5.880 6.40e-09 ***
## pct_esp_vert_high_2011b250 0.056012   0.008764   6.391 3.03e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.759 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.088,  Adjusted R-squared:  0.08535 
## F-statistic: 33.15 on 2 and 687 DF,  p-value: 1.811e-14

5.1.3.3.3 Buffer 500m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_diff_high_2011.2017b500, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_diff_high_2011.2017b500, na.rm = TRUE),
    sd = sd(pct_esp_vert_diff_high_2011.2017b500, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_diff_high_2011.2017b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1   70.9   70.89   25.03 7.18e-07 ***
## Residuals            688 1948.5    2.83                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b500 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b500 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7769 -1.1369 -0.3169  0.7417  7.8431 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               2.80691    0.07665  36.618  < 2e-16 ***
## gentrified_2016_2011TRUE  0.69848    0.13961   5.003 7.18e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.683 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.0351, Adjusted R-squared:  0.0337 
## F-statistic: 25.03 on 1 and 688 DF,  p-value: 7.176e-07

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b500 ~ gentrified_2016_2011 + pct_esp_vert_high_2011b500, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b500 ~ gentrified_2016_2011 + 
##     pct_esp_vert_high_2011b500, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.2608 -1.0193 -0.2532  0.6779  7.6584 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.793392   0.172767  10.380  < 2e-16 ***
## gentrified_2016_2011TRUE   0.816183   0.136808   5.966 3.90e-09 ***
## pct_esp_vert_high_2011b500 0.056440   0.008682   6.501 1.53e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.635 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.09102,    Adjusted R-squared:  0.08838 
## F-statistic:  34.4 on 2 and 687 DF,  p-value: 5.795e-15

5.1.3.3.4 Buffer 750m

ggplot(drop_na(units::drop_units(bei_df_aoi), gentrified_2016_2011), aes(y=pct_esp_vert_diff_high_2011.2017b750, x=gentrified_2016_2011)) +
  geom_boxplot()

bei_df_aoi %>% 
  units::drop_units() %>%
  drop_na() %>%
  group_by(gentrified_2016_2011) %>%
  summarise(
    count = n(),
    mean = mean(pct_esp_vert_diff_high_2011.2017b750, na.rm = TRUE),
    sd = sd(pct_esp_vert_diff_high_2011.2017b750, na.rm = TRUE)
  )

# Compute the analysis of variance
res.aov <- aov(pct_esp_vert_diff_high_2011.2017b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
# Summary of the analysis
summary(res.aov)

##                       Df Sum Sq Mean Sq F value   Pr(>F)    
## gentrified_2016_2011   1   68.7   68.74   28.66 1.17e-07 ***
## Residuals            688 1649.8    2.40                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 15 observations deleted due to missingness

# Linear model
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b750 ~ gentrified_2016_2011, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b750 ~ gentrified_2016_2011, 
##     data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7009 -1.0767 -0.2509  0.6966  6.7491 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               2.74089    0.07053  38.859  < 2e-16 ***
## gentrified_2016_2011TRUE  0.68781    0.12847   5.354 1.17e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.549 on 688 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:   0.04,  Adjusted R-squared:  0.0386 
## F-statistic: 28.67 on 1 and 688 DF,  p-value: 1.174e-07

# Accounting for UC in 2011
res.lm <- lm(pct_esp_vert_diff_high_2011.2017b750 ~ gentrified_2016_2011 + pct_esp_vert_high_2011b750, data = units::drop_units(bei_df_aoi))
summary(res.lm)

## 
## Call:
## lm(formula = pct_esp_vert_diff_high_2011.2017b750 ~ gentrified_2016_2011 + 
##     pct_esp_vert_high_2011b750, data = units::drop_units(bei_df_aoi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1621 -0.9476 -0.2480  0.6476  6.5868 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.708692   0.164555  10.384  < 2e-16 ***
## gentrified_2016_2011TRUE   0.792147   0.125249   6.325 4.57e-10 ***
## pct_esp_vert_high_2011b750 0.057624   0.008359   6.894 1.23e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.499 on 687 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.1021, Adjusted R-squared:  0.0995 
## F-statistic: 39.06 on 2 and 687 DF,  p-value: < 2.2e-16

6 Summary of association between SES and Built Environment

Below we summarize the trends linking SES and BE (looking at association at the 500m buffer scale and INTERACT study area).

SES metric	BE metric	UC 2011 trend	UI 2011 -> 2016 trend
Pampalon / MAT 2011	Bike lane	-***	x (-** when controlling for UC)
Pampalon / MAT 2011	Greenness	-***	x (-* when controlling for UC)
Pampalon / MAT 2011	Tree canopy	-***	-***
Visible Minority 2011	Bike lane	-***	-***
Visible Minority 2011	Greenness	-***	-* (-** when controlling for UC)
Visible Minority 2011	Tree canopy	-***	-***
Gentrified 2011-2016	Bike lane	x	+***
Gentrified 2011-2016	Greenness	-***	+** (x when controlling for UC)
Gentrified 2011-2016	Tree canopy	-***	+***

x: no significant association
-: negative association (higher SES metric => lower BE metric)
+: positive association (higher SES metric => higher BE metric)

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Controlling for UC in UI models does not change the association trend; even better, non significant association may become (slightly) significant.

7 Exploratory stuff

A few ref to properly plot interactions: - https://stats.oarc.ucla.edu/r/seminars/interactions-r/#s4d - https://cran.r-project.org/web/packages/emmeans/vignettes/sophisticated.html

7.1 Urban conditions and SES

Tackling part 2 of objective #1: a multilevel model to test the reduction in socio-economic inequalities in urban conditions \[𝔼(UC_{ij} \mid X_{ij}) = \beta_{0j} + \beta_{1j} ∗ SES + \beta_{2j} ∗ Time + \beta_{3j} ∗ SES ∗ Time + \epsilon_{ij}\]

Data for bike lanes is available for 2006, 2011 and 2016 whereas greenness (and trees) is only available for 2011 and 2016.

# Get SCOREMAT
.bei_df_mlm.1 <- bei_df %>%
  select(CT_UID, starts_with("wSCOREMAT")) %>%
  pivot_longer(cols = starts_with("wSCOREMAT"), 
               names_to = "Year", names_prefix = "wSCOREMAT.", 
               values_to = "wSCOREMAT", values_drop_na = TRUE)
# Get vis_minority
.bei_df_mlm.2 <- bei_df %>%
  select(CT_UID, starts_with("vis_minority")) %>%
  pivot_longer(cols = starts_with("vis_minority"), 
               names_to = "Year", names_prefix = "vis_minority_", 
               values_to = "vis_minority", values_drop_na = TRUE)
# Get gentrified flag
.bei_df_mlm.3 <- bei_df %>%
  select(CT_UID, starts_with("gentrified_")) %>%
  pivot_longer(cols = starts_with("gentrified_"), 
               names_to = "Year_span", names_prefix = "gentrified_", 
               values_to = "gentrified", values_drop_na = TRUE) %>%
  extract(Year_span, "Year", "(\\d{4})_*")
# Get bike_lane
.bei_df_mlm.4 <- bei_df %>%
  select(CT_UID, starts_with("Bike_lane.by.street")) %>%
  pivot_longer(cols = starts_with("Bike_lane.by.street"), 
               names_to = "Var", names_prefix = "Bike_lane.by.street.",
               values_to = "bike_lane.by.street", values_drop_na = TRUE) %>%
  extract("Var", c("Year", "Spatial.scale"), "(\\d{4})([[:alnum:]]+)")
# Get Canopy
.bei_df_mlm.5 <- bei_df %>%
  select(CT_UID, starts_with("pct_esp_vert_20")) %>%
  pivot_longer(cols = starts_with("pct_esp_vert_"), 
               names_to = "Var", names_prefix = "pct_esp_vert_",
               values_to = "pct_esp_vert", values_drop_na = TRUE) %>%
  extract("Var", c("Year", "Spatial.scale"), "(\\d{4})([[:alnum:]]+)") %>%
  filter(Year %in% c('2011', '2017')) %>%
  mutate(Year = case_when(Year == '2017' ~ '2016',
                          TRUE ~ Year))
# Get Canopy
.bei_df_mlm.6 <- bei_df %>%
  select(CT_UID, starts_with("pct_esp_vert_high_20")) %>%
  pivot_longer(cols = starts_with("pct_esp_vert_high"), 
               names_to = "Var", names_prefix = "pct_esp_vert_high_",
               values_to = "pct_esp_vert_high", values_drop_na = TRUE) %>%
  extract("Var", c("Year", "Spatial.scale"), "(\\d{4})([[:alnum:]]+)") %>%
  filter(Year %in% c('2011', '2017')) %>%
  mutate(Year = case_when(Year == '2017' ~ '2016',
                          TRUE ~ Year))

# Combine all subsets
bei_df_mlm <- bei_df %>%
  select(CT_UID, interact_aoi, Population) %>%
  full_join(.bei_df_mlm.1, by=c("CT_UID")) %>%
  full_join(.bei_df_mlm.2, by=c("CT_UID", "Year")) %>%
  full_join(.bei_df_mlm.3, by=c("CT_UID", "Year")) %>%
  full_join(.bei_df_mlm.4, by=c("CT_UID", "Year")) %>%
  full_join(.bei_df_mlm.5, by=c("CT_UID", "Year", "Spatial.scale")) %>%
  full_join(.bei_df_mlm.6, by=c("CT_UID", "Year", "Spatial.scale")) %>%
  units::drop_units() %>%
  #mutate(Year = factor(Year, levels = c("2006", "2011", "2016")))
  mutate(Year = as.numeric(Year))

7.1.1 Model 1: SCOREMAT & bike lanes

Here, SES is measured through Pampalon Material Score and urban conditions as the ratio of bike lanes to streets within 500m buffers around CT for the INTERACT study area. Using 2006 as reference year.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi) %>%
  mutate(delta.Year = Year - 2006)

res.mlm.1 <- lmer(bike_lane.by.street ~ wSCOREMAT*delta.Year + (1 | CT_UID), data = mlm_data, REML = FALSE)
summary(res.mlm.1)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bike_lane.by.street ~ wSCOREMAT * delta.Year + (1 | CT_UID)
##    Data: mlm_data
## 
##      AIC      BIC   logLik deviance df.resid 
##  12079.8  12113.6  -6033.9  12067.8     2066 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6767 -0.5590  0.0050  0.5021  3.9861 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 21.75    4.664   
##  Residual             10.14    3.185   
## Number of obs: 2072, groups:  CT_UID, 693
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)            6.12580    0.20903  29.306
## wSCOREMAT            -29.66509    4.78235  -6.203
## delta.Year             0.60638    0.01715  35.363
## wSCOREMAT:delta.Year  -1.93134    0.44834  -4.308
## 
## Correlation of Fixed Effects:
##             (Intr) wSCOREMAT dlt.Yr
## wSCOREMAT   -0.029                 
## delta.Year  -0.410 -0.002          
## wSCOREMAT:.  0.004 -0.506    -0.007

# Plot interactions
plot_model(res.mlm.1, type = 'int')

emmip(res.mlm.1, delta.Year~wSCOREMAT, at = list(delta.Year = c(0, 5, 10), wSCOREMAT=seq(-.2, .2, by=.01)), CI=TRUE)

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

7.1.2 Model 2: Visibility minority & bike lanes

Here, SES is measured through the % of visible minority population and urban conditions as the ratio of bike lanes to streets within 500m buffers around CT for the INTERACT study area.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi) %>%
  mutate(delta.Year = Year - 2006)

res.mlm.2 <- lmer(bike_lane.by.street ~ vis_minority*delta.Year + (1 | CT_UID), data = mlm_data, REML = FALSE)
summary(res.mlm.2)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bike_lane.by.street ~ vis_minority * delta.Year + (1 | CT_UID)
##    Data: mlm_data
## 
##      AIC      BIC   logLik deviance df.resid 
##  11934.5  11968.3  -5961.2  11922.5     2059 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6629 -0.5617 -0.0067  0.4941  3.7786 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 19.187   4.380   
##  Residual              9.933   3.152   
## Number of obs: 2065, groups:  CT_UID, 691
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              8.270127   0.329009  25.137
## vis_minority            -0.106487   0.012148  -8.766
## delta.Year               0.845371   0.032637  25.902
## vis_minority:delta.Year -0.005666   0.001114  -5.086
## 
## Correlation of Fixed Effects:
##             (Intr) vs_mnr dlt.Yr
## vis_minorty -0.795              
## delta.Year  -0.340  0.245       
## vs_mnrty:.Y  0.437 -0.537 -0.821

# Plot interactions
plot_model(res.mlm.2, type = 'int')

7.1.3 Model 3: SCOREMAT, visible minority & bike lanes

Here, SES is measured through Pampalon Material Score combined with the % of visible minority population and urban conditions as the ratio of bike lanes to streets within 500m buffers around CT for the INTERACT study area.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi) %>%
  mutate(delta.Year = Year - 2006)

res.mlm.3 <- lmer(bike_lane.by.street ~ wSCOREMAT*delta.Year + vis_minority*delta.Year + (1 | CT_UID), data = mlm_data, REML = FALSE)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

summary(res.mlm.3)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bike_lane.by.street ~ wSCOREMAT * delta.Year + vis_minority *  
##     delta.Year + (1 | CT_UID)
##    Data: mlm_data
## 
##      AIC      BIC   logLik deviance df.resid 
##  11920.3  11965.4  -5952.1  11904.3     2057 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6895 -0.5686  0.0009  0.4971  3.8313 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 19.285   4.391   
##  Residual              9.792   3.129   
## Number of obs: 2065, groups:  CT_UID, 691
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)               7.877530   0.351323  22.422
## wSCOREMAT               -16.282727   5.178330  -3.144
## delta.Year                0.812943   0.036822  22.077
## vis_minority             -0.086567   0.013567  -6.381
## wSCOREMAT:delta.Year     -0.491004   0.518795  -0.946
## delta.Year:vis_minority  -0.005020   0.001297  -3.871
## 
## Correlation of Fixed Effects:
##             (Intr) wSCOREMAT dlt.Yr vs_mnr wSCOREMAT:
## wSCOREMAT    0.355                                   
## delta.Year  -0.324 -0.139                            
## vis_minorty -0.823 -0.449     0.242                  
## wSCOREMAT:. -0.174 -0.499     0.463  0.203           
## dlt.Yr:vs_m  0.438  0.260    -0.858 -0.514 -0.522    
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

# Plot interactions
plot_model(res.mlm.3, type = 'pred', terms = c('wSCOREMAT', 'vis_minority', 'delta.Year'))

7.1.4 Model 4: SCOREMAT, visible minority & esp_vert

Here, SES is measured through Pampalon Material Score combined to % of visible minority population and urban conditions as the % of greenness within 500m buffers around CT for the INTERACT study area.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi) %>%
  mutate(delta.Year = Year - 2006)

res.mlm.4 <- lmer(pct_esp_vert ~ wSCOREMAT*delta.Year + vis_minority*delta.Year + (1 | CT_UID), data = mlm_data)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

summary(res.mlm.4)

## Linear mixed model fit by REML ['lmerMod']
## Formula: pct_esp_vert ~ wSCOREMAT * delta.Year + vis_minority * delta.Year +  
##     (1 | CT_UID)
##    Data: mlm_data
## 
## REML criterion at convergence: 8727.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.1003 -0.3850  0.0039  0.3822  8.3635 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 143.15   11.965  
##  Residual               3.66    1.913  
## Number of obs: 1378, groups:  CT_UID, 691
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)              32.779289   0.693038  47.298
## wSCOREMAT               -36.172408   7.773713  -4.653
## delta.Year                0.845882   0.049396  17.125
## vis_minority             -0.031168   0.021497  -1.450
## wSCOREMAT:delta.Year      0.003342   0.626849   0.005
## delta.Year:vis_minority  -0.002711   0.001556  -1.743
## 
## Correlation of Fixed Effects:
##             (Intr) wSCOREMAT dlt.Yr vs_mnr wSCOREMAT:
## wSCOREMAT    0.156                                   
## delta.Year  -0.297 -0.026                            
## vis_minorty -0.707 -0.325     0.187                  
## wSCOREMAT:. -0.162 -0.491     0.477  0.196           
## dlt.Yr:vs_m  0.401  0.232    -0.861 -0.482 -0.552    
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

# Plot interactions
plot_model(res.mlm.4, type = 'pred', terms = c('wSCOREMAT', 'vis_minority', 'delta.Year'))

7.1.5 Model 5: SCOREMAT, vis_minority & esp_vert_high

Here, SES is measured through Pampalon Material Score combined with % of visible minority population and urban conditions as the % of trees within 500m buffers around CT for the INTERACT study area.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi) %>%
  mutate(delta.Year = Year - 2006)

res.mlm.5 <- lmer(pct_esp_vert_high ~ wSCOREMAT*delta.Year + vis_minority*delta.Year + (1 | CT_UID), data = mlm_data)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

summary(res.mlm.5)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## pct_esp_vert_high ~ wSCOREMAT * delta.Year + vis_minority * delta.Year +  
##     (1 | CT_UID)
##    Data: mlm_data
## 
## REML criterion at convergence: 7254.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.98339 -0.41949  0.00332  0.40454  3.06727 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 46.920   6.850   
##  Residual              1.302   1.141   
## Number of obs: 1378, groups:  CT_UID, 691
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)              1.540e+01  4.042e-01  38.100
## wSCOREMAT               -1.923e+01  4.604e+00  -4.177
## delta.Year               6.103e-01  2.942e-02  20.744
## vis_minority            -3.527e-02  1.269e-02  -2.779
## wSCOREMAT:delta.Year    -2.492e+00  3.738e-01  -6.667
## delta.Year:vis_minority -5.059e-04  9.277e-04  -0.545
## 
## Correlation of Fixed Effects:
##             (Intr) wSCOREMAT dlt.Yr vs_mnr wSCOREMAT:
## wSCOREMAT    0.162                                   
## delta.Year  -0.308 -0.030                            
## vis_minorty -0.716 -0.331     0.195                  
## wSCOREMAT:. -0.167 -0.496     0.477  0.200           
## dlt.Yr:vs_m  0.411  0.235    -0.862 -0.487 -0.551    
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

# Plot interactions
plot_model(res.mlm.5, type = 'pred', terms = c('wSCOREMAT', 'vis_minority', 'delta.Year'))

7.1.6 Model 6: SCOREMAT & bike lanes, random slope

Here, SES is measured through Pampalon Material Score and urban conditions as the ratio of bike lanes to streets within 500m buffers around CT for the INTERACT study area. Modeled as random slope, not just random intercept.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi)  %>%
  mutate(delta.Year = Year - 2006)

res.mlm.6 <- lmer(bike_lane.by.street ~ wSCOREMAT*delta.Year + (delta.Year | CT_UID), data = mlm_data, REML = FALSE)
summary(res.mlm.6)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bike_lane.by.street ~ wSCOREMAT * delta.Year + (delta.Year |  
##     CT_UID)
##    Data: mlm_data
## 
##      AIC      BIC   logLik deviance df.resid 
##  11900.1  11945.2  -5942.1  11884.1     2064 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.00916 -0.40756  0.00695  0.36142  2.84899 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  CT_UID   (Intercept) 20.2983  4.5054        
##           delta.Year   0.1717  0.4144   -0.08
##  Residual              5.8716  2.4231        
## Number of obs: 2072, groups:  CT_UID, 693
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)            6.12408    0.19092  32.076
## wSCOREMAT            -24.50034    4.37362  -5.602
## delta.Year             0.60615    0.02046  29.621
## wSCOREMAT:delta.Year  -2.18277    0.51916  -4.204
## 
## Correlation of Fixed Effects:
##             (Intr) wSCOREMAT dlt.Yr
## wSCOREMAT   -0.029                 
## delta.Year  -0.272 -0.001          
## wSCOREMAT:. -0.004 -0.418    -0.007

# Check if model with random slope significantly improves over random intercept (see http://www.bodowinter.com/tutorial/bw_LME_tutorial1.pdf)
res.mlm.6.null <- lmer(bike_lane.by.street ~ wSCOREMAT*delta.Year + (1 | CT_UID), data = mlm_data, REML = FALSE)
anova(res.mlm.6.null, res.mlm.6)

7.1.7 Model 7: SCOREMAT, vis_minority & bike lanes, three-way iteraction

Here, SES is measured through Pampalon Material and % of visible minority and urban conditions as the ratio of bike lanes to streets within 500m buffers around CT for the INTERACT study area. Here, the three-way interaction is modeled, contrary to model 3, where the same variables are considered through two-way interaction (Year & SCOREMAT and Year & vis_minority).

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi)  %>%
  mutate(delta.Year = Year - 2006)

res.mlm.7 <- lmer(bike_lane.by.street ~ wSCOREMAT*vis_minority*delta.Year + (1 | CT_UID), data = mlm_data, REML = FALSE)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

summary(res.mlm.7)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bike_lane.by.street ~ wSCOREMAT * vis_minority * delta.Year +  
##     (1 | CT_UID)
##    Data: mlm_data
## 
##      AIC      BIC   logLik deviance df.resid 
##  11923.4  11979.8  -5951.7  11903.4     2055 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7095 -0.5622  0.0039  0.4893  3.8351 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 19.305   4.394   
##  Residual              9.782   3.128   
## Number of obs: 2065, groups:  CT_UID, 691
## 
## Fixed effects:
##                                     Estimate Std. Error t value
## (Intercept)                         7.921273   0.356795  22.201
## wSCOREMAT                         -20.930417   7.600014  -2.754
## vis_minority                       -0.091188   0.014800  -6.161
## delta.Year                          0.811819   0.037344  21.739
## wSCOREMAT:vis_minority              0.202550   0.241270   0.840
## wSCOREMAT:delta.Year               -0.519784   0.818174  -0.635
## vis_minority:delta.Year            -0.004881   0.001427  -3.421
## wSCOREMAT:vis_minority:delta.Year  -0.003483   0.024586  -0.142
## 
## Correlation of Fixed Effects:
##              (Intr) wSCOREMAT vs_mnr dlt.Yr wSCOREMAT:v_ wSCOREMAT:. vs_:.Y
## wSCOREMAT     0.115                                                        
## vis_minorty  -0.811  0.010                                                 
## delta.Year   -0.337 -0.019     0.264                                       
## wSCOREMAT:v_  0.168 -0.732    -0.396 -0.099                                
## wSCOREMAT:.  -0.026 -0.435    -0.034  0.162  0.288                         
## vs_mnrty:.Y   0.451 -0.037    -0.550 -0.839  0.265        0.008            
## wSCOREMAT:_: -0.132  0.414     0.263  0.169 -0.552       -0.758      -0.416
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

# Plot interactions
plot_model(res.mlm.7, type = 'pred', terms = c('wSCOREMAT', 'vis_minority', 'delta.Year'))

# Check if model with three-way interaction significantly improves over two-way interactions (model 3)
anova(res.mlm.3, res.mlm.7)

7.1.8 Model 8: Gentrification & bike lanes

Here, SES is measured through gentrification status and urban conditions as the ratio of bike lanes to streets within 500m buffers around CT for the INTERACT study area. Using 2006 as reference year.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi)  %>%
  mutate(delta.Year = Year - 2006)

res.mlm.8 <- lmer(bike_lane.by.street ~ factor(gentrified)*delta.Year + (1 | CT_UID), data = mlm_data, REML = FALSE)
summary(res.mlm.8)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bike_lane.by.street ~ factor(gentrified) * delta.Year + (1 |  
##     CT_UID)
##    Data: mlm_data
## 
##      AIC      BIC   logLik deviance df.resid 
##  11673.1  11706.7  -5830.6  11661.1     1991 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8787 -0.5363 -0.0127  0.4797  3.9410 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 23.081   4.804   
##  Residual              9.962   3.156   
## Number of obs: 1997, groups:  CT_UID, 691
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                        6.74580    0.23623  28.556
## factor(gentrified)TRUE            -2.40404    0.29387  -8.181
## delta.Year                         0.47163    0.02195  21.485
## factor(gentrified)TRUE:delta.Year  0.44069    0.04161  10.590
## 
## Correlation of Fixed Effects:
##             (Intr) fc()TRUE dlt.Yr
## fctr(g)TRUE -0.399                
## delta.Year  -0.487  0.432         
## fc()TRUE:.Y  0.291 -0.710   -0.600

# Plot interactions
plot_model(res.mlm.8, type = 'pred', terms = c('gentrified', 'delta.Year'))

7.1.9 Model 9: Gentrification & esp_vert

Here, SES is measured through gentrification status and urban conditions as the ratio of bike lanes to streets within 500m buffers around CT for the INTERACT study area. Using 2006 as reference year.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi)  %>%
  mutate(delta.Year = Year - 2006)

res.mlm.9 <- lmer(pct_esp_vert ~ factor(gentrified)*delta.Year + (1 | CT_UID), data = mlm_data, REML = FALSE)
summary(res.mlm.9)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: pct_esp_vert ~ factor(gentrified) * delta.Year + (1 | CT_UID)
##    Data: mlm_data
## 
##      AIC      BIC   logLik deviance df.resid 
##   8691.4   8722.7  -4339.7   8679.4     1357 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.2575 -0.3735 -0.0241  0.3650  8.5444 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 149.043  12.208  
##  Residual               3.732   1.932  
## Number of obs: 1363, groups:  CT_UID, 689
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                       32.44059    0.51570  62.906
## factor(gentrified)TRUE            -1.89689    0.46658  -4.066
## delta.Year                         0.71055    0.02757  25.769
## factor(gentrified)TRUE:delta.Year  0.21239    0.05747   3.696
## 
## Correlation of Fixed Effects:
##             (Intr) fc()TRUE dlt.Yr
## fctr(g)TRUE -0.285                
## delta.Year  -0.405  0.600         
## fc()TRUE:.Y  0.261 -0.920   -0.646

# Plot interactions
plot_model(res.mlm.9, type = 'pred', terms = c('gentrified', 'delta.Year'))

7.1.10 Model 10: Gentrification & esp_vert

Here, SES is measured through gentrification status and urban conditions as the ratio of bike lanes to streets within 500m buffers around CT for the INTERACT study area. Using 2006 as reference year.

# Define multilevel model | buf 500m / INTERACT study area
mlm_data <- bei_df_mlm %>%
  filter(Spatial.scale == 'b500' & interact_aoi)  %>%
  mutate(delta.Year = Year - 2006)

res.mlm.10 <- lmer(pct_esp_vert_high ~ factor(gentrified)*delta.Year + (1 | CT_UID), data = mlm_data, REML = FALSE)
summary(res.mlm.10)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: pct_esp_vert_high ~ factor(gentrified) * delta.Year + (1 | CT_UID)
##    Data: mlm_data
## 
##      AIC      BIC   logLik deviance df.resid 
##   7330.6   7361.9  -3659.3   7318.6     1357 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.97493 -0.40607 -0.00555  0.40257  3.14651 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CT_UID   (Intercept) 54.531   7.385   
##  Residual              1.385   1.177   
## Number of obs: 1363, groups:  CT_UID, 689
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                       14.96985    0.31234  47.927
## factor(gentrified)TRUE            -2.08819    0.28420  -7.348
## delta.Year                         0.52563    0.01680  31.294
## factor(gentrified)TRUE:delta.Year  0.25392    0.03501   7.253
## 
## Correlation of Fixed Effects:
##             (Intr) fc()TRUE dlt.Yr
## fctr(g)TRUE -0.287                
## delta.Year  -0.408  0.600         
## fc()TRUE:.Y  0.262 -0.920   -0.646

# Plot interactions
plot_model(res.mlm.10, type = 'pred', terms = c('gentrified', 'delta.Year'))